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The  objective  of  this  thesis  is  to  develop  an 
interactive  computer  progr6un  that  will  analyze  and 
synthesize  radiative  heat  transfer  in  longitudinal  fins. 

The  analysis  procedure  determines  the  aunount  of  heat 
transferred  from  the  fin  given  the  fin  base  temperature,  fin 
dimensions  and  thermal  properties.  The  synthesis  procedure 
is  the  converse  problem:  it  determines  the  size  of  the  fin 
required  to  dissipate  a  specified  amount  of  heat  given  the 
thermal  characteristics  of  the  fin.  In  addition,  the 
program  is  capable  of  performing  the  analysis/synthesis  of 
three  fin  profiles  (rectangular,  trapezoidal  and  triangular) 
in  two  environments  (free  space  and  non-free  space).  Free 
space  is  considered  as  the  absence  of  external  heat  sources 
or  interception  of  the  heat  dissipated  by  the  fin  whereas 
non-free  space  includes  the  effect  of  external  heat  sources 
and  neighboring  structures.  A  theoretical  analysis  of  heat 
transfer  from  radiating  longitudinal  fins  will  be  presented 
along  with  a  user  oriented  computer  program.  Finally, 
detailed  examples  will  be  provided  to  illustrate  the 
different  types  of  problems,  profiles  and  environments. 
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I .  INTRODUCTION 


A.  DEFINITIONS 

The  American  Heritage  Dictionary  of  the  English  Language 
defines  heat  as  "  a  form  of  energy  associated  with  the 
motion  of  atoms  or  molecules  in  solids  and  capable  of  being 
transmitted  through  solid  and  fluid  media  by  conduction, 
through  fluid  media  by  convection,  and  through  empty  space 
by  radiation  "  [Ref.  1:  p.  608],  The  transfer  of  heat 
within  a  homogeneous  substance  because  of  a  temperature 
differential  is  called  conduction.  Convection  is  the 
transfer  of  heat  between  a  solid  body  and  a  contiguous, 
moving  fluid  at  a  different  temperature.  Radiation  is  the 
transfer  of  heat  by  electromagnetic  waves  in  a  medium  or  a 
vacuum. 

The  American  Heritage  Dictionary  of  the  English  Language 
defines  a  fin  as  "a  projecting  vane  used  for  cooling,  as  on 
a  radiator  or  engine  cylinder"  [Ref.  1:  p.  492].  The  use  of 
fins  is  a  commonly  used  technique  to  increase  exposed 
surface  area  and  hence  increase  the  rate  of  heat  transfer 
from  a  body  to  its  surroundings. 

B.  EXAMPLES  OF  HEAT  TRANSFER  BY  FINS 

Applications  of  heat  transfer  augmentation  through  the 
use  of  fins  range  from  motorcycle  engine  blocks  to  space 
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stations.  A  motorcycle's  internal  combustion  gasoline 
engine  generates  a  large  amount  of  heat  that  must  be 
dissipated  in  order  for  the  engine  to  remain  within  the 
permissible  temperature  limits  of  the  cylinder  liners.  The 
engine  casing  does  not  provide  sufficient  surface  area  to 
dissipate  this  heat;  therefore,  external  fins  are  installed 
to  increase  the  heat  transfer  as  indicated  in  Figures  lA  and 
IB  [Ref.  2:  pp.  1,4,56].  Similarly  the  Space  Station, 
"Freedom",  will  require  hundreds  of  square  meters  of  surface 
area  to  dissipate  the  heat  generated  by  its  electronic 
equipment  as  shown  in  Figure  2  [Ref.  3:  p.  29]. 


Figure  lA  Fins  on  Motorcycle  Engine 


Figure  IB  Fins  on  Motorcycle  Engine 


Figure  2  Radiators  on  Space  Station 

C.  PROFILES 

Fins  are  fabricated  in  many  sizes  and  shapes  [Ref. 

4:  pp.  102  -  103],  [Ref.  5:  pp.  3-22];  however,  for  this 
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thesis  only  three  common  shapes  of  longitudinal  fins  are 
investigated.  These  are  the  fins  of  rectangular. 
trapezoidal  and  triangular  profile.  Longitudinal  in  the 
foregoing  sense  means  "placed  or  running  lengthwise"  [Ref. 

1;  p.  769];  that  is  to  say,  the  analysis  is  simplified  to  a 
one  dimensional  variation  of  temperature  defined  along  the 
"x"  axis  and  a  uniform  distribution  of  temperature  along  the 
other  two  axes.  Figure  3  shows  these  profiles. 


Figure  3  Longitudinal  Fin  Profiles 
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D.  ENVIRONMENTS 

The  three  types  of  fins  will  be  analyzed  with  respect  to 
two  environments:  free  space  and  non-free  space.  Free 
space  in  this  sense  means  free  of  external  heat  sources  as 
seen  by  the  fin  or  radiator  or  free  of  nearby  surfaces  that 
intercept  radiating  energy  from  the  fin  .  Non-free  space 
includes  the  effect  of  external  heat  sources.  For  ex6anple/ 
a  fin  exposed  to  sunlight  will  have  an  external  solar  heat 
flux  incident  upon  it  in  addition  to  the  heat  conducted 
through  the  base  of  the  fin.  A  fin  in  the  shade  will  have 
no  such  external  source  and  therefore  may  radiate  into  "free 
space".  For  the  non-free  space  environment,  the  term 
"external  heat"  will  also  include  the  following: 

1.  Radiative  heat  incident  on  the  fin  from  external  heat 
sources . 

2.  Radiative  heat  from  the  fin  intercepted  by  nearby 
surfaces  such  as  solar  panels  or  other  dissipating 
fins  or  the  vehicle  structure  itself. 

The  radiative  transfer  between  a  fin  and  its  environment 
can  be  quite  complex  depending  on  the  number  of  heat 
sources,  the  geometric  shapes  of  the  heat  sources  and  their 
special  arrangements  (view  factors)  [Ref.  6:  pp.  1  -  7]. 

An  idea  of  the  free  space  and  non-free  space 
environments  is  shown  in  Figure  4 . 

E .  PROBLEM  TYPES 

Two  types  of  problems  will  be  investigated.  Analysis  is 
the  problem  of  determining  the  eunount  of  heat  transferred  to 
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Figure  4  Free  Space  Vs  Non-Free  Space  Environments 


the  environment  from  a  fin  with  given  dimensions,  base 
temperature,  external  heat  sources  nearby  and  particular 
thermal  properties.  Synthesis  is  the  converse  problem; 
given  a  specified  amount  of  heat  to  be  dissipated  at  a 
certain  base  temperature  and  with  certain  thermal 
properties,  the  fin  is  sized  so  as  to  efficiently  transmit 
this  heat  to  the  environment.  These  input /output 
relationships  are  shown  in  Figure  5. 


F.  STEADY-STATE  VS  NON-STEADT-STATE 

Only  a  steady-state  is  considered.  Within  this  static 
context,  the  temperature  does  not  change  with  time. 
Mathematically  this  means  that  any  derivatives  of 
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BASE  TEMPERATURE  - 

FIN  DIMENSIONS 

EDGE  TEMPERATURE 

CONDUCTIVITY 

— ►  HEAT  DISSIPATED 

ABSORPTIVITY 

FIN  EFFICIENCY 

EMISSIVITY 

EXTERNAL  HEAT  INPUT  - 

EFFECT  OF  NEARBY  BODIES 

ANALYSIS 

HEAT  DISSIPATION  REQUIRED  - 

BASE  TEMPERATURE 

CONDUCTIVITY 

EDGE  TEMPERATURE 

ABSORPTIVITY 

— ►  FIN  DIMENSIONS 

EMISSIVITY 

FIN  EFFICIENCY 

EXTERNAL  HEAT  INPUT  - 

EFFECT  OF  NEARBY  BODIES 

SYNTHESIS 

Figure  5  Analysis  vs  Synthesis  Problems 


temperature  with  respect  to  time  are  zero.  This  is  in 
contrast  to  the  non- steady- state  or  dynamic  situation  where 
temperature  is  a  function  of  time.  A  satellite  emerging 
from  the  Earth's  shadow  would  experience  a  warming  period  of 
increasing  temperature  (non-steady-state)  until  it  reaches 
thermal  equilibrium  with  respect  to  the  space  environment  at 
which  time,  the  temperature  would  remain  constant 
(steady-state) . 
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G.  OVERVIEW 


In  summary,  this  thesis  will  analyze  longitudinal  fins 
of  three  basic  profiles  (rectangular,  trapezoidal  and 
triangular)  in  two  environments  (free  space  and  non-free 
space)  to  obtain  the  solution  to  two  types  of  problems 
(analysis  and  synthesis).  After  a  theoretical  analysis  is 
developed,  a  computer  progreun  will  be  presented  that  will 
treat  these  situations  through  an  interactive  series  of  user 
oriented  menus.  Finally,  detailed  examples  of  typical 
applications  will  be  presented.  An  overview  is  shown  in 
Figure  6. 


FREE  SPACE— I 
NON-FREE  space-4 


RECTANGLE- 
TRAPEZOID — 
TRIANGLE- 

ANALYSIS- 

SYNTHESIS—* 


THEORETICAL 

ANALYSIS 


ALGORITHM  -fr 


COMPUTER  - EXAMPLES 

PROGRAM 


Figure  6  Overview 
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II.  ANALYSIS  OF  RADIATIVE  HEAT  TRANSFER  IN  LONGITUDINAL 

FINS 


A.  PRELIMINARY 

1 .  Restrictions 

The  analysis  will  be  restricted  to  an  environment 
where  the  absence  of  a  fluid  medium  (vacuum)  precludes  heat 
transfer  by  convection.  Here,  only  the  non-free  case  will 
be  analyzed  for  each  of  the  profiles  (rectangular, 
trapezoidal,  triangular)  and  the  two  types  problems 
(analysis,  synthesis)  since  it  is  the  more  general  case  and 
includes  the  free  space  case. 

2.  General  Heat  Balance  Equation 

The  general  heat  balance  equation  as  derived  from 
the  principle  of  the  conservation  of  energy  states  that  the 
amount  of  heat  gained  by  an  object  from  its  surroundings  is 
equal  to  the  heat  lost  plus  any  heat  stored  by  the  object. 

HEAT  GAINED  =  HEAT  LOST  +  HEAT  STORED  [1] 

3.  Heat  Transfer  by  Conduction 

The  transfer  of  heat  within  a  homogeneous  substance 
because  of  a  temperature  differential  is  called  conduction. 
The  simplest  case  of  one  dimensional  heat  flow  by  conduction 
in  a  solid  is  shown  in  Figure  7. 
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where : 

Q  =  heat  flow  (W) 

k  =  thermal  conductivity  which  is  a  property  of  the 
material  (W/m-°K) 

X  =  length  coordinate  over  which  the  energy  is 
transferred  (m) 

A  =  area  normal  to  the  direction  of  heat  flow  (m^) 
dx  =  differential  element  along  x  (m) 

dT  =  the  temperature  gradient  in  the  x  direction  (‘^K/m) 
dx 

Tg  =  temperature  of  "hot"  end  (°K) 

Tq  =  temperature  of  "cold"  end  (°K) 

The  transfer  of  energy  by  conduction  is  described  by 
Fourier's  law  [Ref.  7;  p.  26]: 


dT  [2] 

Q  =  -  k  A  — 
dx 

Equation  [2]  states  that  the  amount  of  energy 
transferred  is  proportional  to  the  temperature  differential 
along  x  and  the  cross-sectional  area  normal  to  the  direction 


10 


of  heat  flow.  The  constant  of  proportionality,  k,  is  called 
the  thermal  conductivity  and  depends  on  the  material.  The 
minus  sign  conforms  to  the  convention  that  heat  flows  from  a 
hot  source  to  a  cold  sink  and  assures  a  positive  heat  flow 
in  the  presence  of  a  negative  temperature  gradient.  More 
complex  equations  involving  two  and  three  dimensions, 
steady-state  and  transient  heat  flow  in  different  coordinate 
systems  (rectangular,  cylindrical  and  spherical)  can  be 
derived  [Ref.  7:  p.  17]  (the  one  dimensional  case  is  a 
considerable  simplification)  but  will  not  be  given  here. 

The  thermal  conductivities  of  some  spacecraft  materials  are 
given  in  Table  I  [Ref.  8:  p.  269]. 

TABLE  I  THERMAL  CONDUCTIVITY 


Material  (at  298  °K) 

k  (W/m-OR) 

Aluminum 

210 

Aluminum  alloys 

117  -  175 

Magnesium 

157 

Magnesium  alloys 

52  -  111 

Titanium 

21 

Stainless  steel 

16.2 

Teflon 

0.25 

4.  Heat  Transfer  by  Radiation 

Radiation  is  the  transfer  of  heat  by  electromagnetic 
waves  in  a  medium  or  a  vacuum.  A  simple  example  is  shown  in 
Figure  8. 


where: 

Q  =  amount  of  heat  radiated  (W) 

e  =  emissivity  of  the  surface  (0  s  e  £l)  (dimensionless) 
a  -  Stef an-Boltzman  constant  (5.6697x10-®  W/m^-K^) 

S  =  surface  area  for  energy  transfer  (m^) 

T  =  temperature  of  the  object  (°K) 

Tg  =  temperature  of  the  environment  (°K) 

The  transfer  of  energy  by  radiation  (to  free  space 
at  T  =  0  °K)  is  described  by  Stef an-Boltzman  Law  [Ref.  9:  p. 
18] : 

Q  =  a  e  S  [3] 

Note  the  presence  of  the  highly  nonlinear  term. 
The  emissivity  of  the  surface  is  determined  by  the  type  of 
material.  Table  II  lists  the  eroissivities  of  some  typical 
spacecraft  materials  [Ref.  8:  p.  275]. 
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TABLE  II  EMISSIVITY 


Material  (9.3;im,310  °K) 

Emissivity 

Black  paint 

0.90 

White  paint 

0.90 

Graphite /epoxy 

0.85 

Launpblack 

0.84 

Solar  cells 

0.82 

Optical  solar  reflector 

0.80 

Anodized  aluminum 

0.80 

Aluminized  kapton 

0.60 

Tiodized  titanium 

0.60 

Aluminum 

0.06 

Gold 

0.03 

B.  RECTANGULAR  FIN 
1 .  Geometry 

The  geometry,  terminology  and  coordinate  system  for 
the  longitudinal  fin  of  rectangular  profile  are  shown  in 
Figure  9: 


Figure  9  Heat  Transfer  in  Rectangular  Fin 
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Here; 


Q  =  amount  of  heat  transfer  (W) 

X  =  height  coordinate  along  the  fin  with  the  positive 
orientation  from  base  to  tip  (m) 

dx  =  differential  element  (m) 

T(x)  =  fin  temperature  at  x  (°K) 

Tq  =  temperature  at  fin  base  (°K) 

Tg  =  temperature  at  fin  tip  (®K) 

Tg  =  temperature  of  environment  (for  free  space, 

Tg  =  0°K) 

L  =  length  of  fin  (m) 
b  =  height  of  fin  (m) 

6  =  thickness  or  width  of  fin  (m) 

A  =  cross-sectional  area  (m^) 

2 .  Assumptions 

The  following  assumptions  were  made  to  simplify  the 
mathematical  development  [Ref.  10:  p.  196]: 

1.  The  temperature  surrounding  the  fin  is  uniform. 

2.  Steady-state. 

3.  No  bond  resistance  exists  between  the  fin  and  the  prime 
surface. 

4.  There  is  no  heat  transfer  through  the  fin  tip. 

6.  The  fin  base  temperature  is  uniform. 

7.  The  fin  width  or  thickness  is  small  compared  with  the 
height  and  length  of  the  fin. 

8.  The  thermal  conductivity  of  the  fin  material  is 
isotropic  and  not  a  function  of  the  temperature. 
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3.  Analysis 

For  each  differential  element  of  the  fin,  the 
difference  between  heat  entering  and  leaving  by  conduction 
is 

d 

dQ  =  k  — 
dx 

the  heat  dissipated  by  radiation  is 

dQR  *  2  a  €  L  t4  dx  [5] 


dT 
A  — 
dx 


dx 


[4] 


and  the  heat  incident  from  external  sources  is 


dQg  =  E  a  L  dx  [6] 

where : 

E  =  external  heat  input  (W/m^) 
a  *  surface  absorptivity  (dimensionless) 

The  absorptivities  of  some  typical  spacecraft  materials 
are  listed  in  Table  III. 

TABLE  III  ABSORPTIVITY 


Material  (9. 3pm, 310  °K) 

Absorptivity 

Black  paint 

0.90 

Graphite /epoxy 

0.84 

Solar  cells 

0.70 

Tiodized  titanium 

0.60 

Aluminized  kapton 

0.35 

Gold 

0.25 

White  paint 

0.20 

Anodized  aluminum 

0.20 

Aluminum 

0.12 

Optical  solar  reflector 

0.08 

Application  of  the  general  energy  balance  [1]  to  the 
differential  element  in  steady-state  results  in 
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dQ  -  dQj^  -  dQg 


[7] 


Substituting  equations  [4]  through  [6]  into  equation  [7] 
yields 


dx 


2  a  e  L  dx  -  E  a  L  dx 


[8] 


With  A=6L,  K]^  =  2oe  and  K2  =  E  a  in  equation  [8], 

simplification  gives 


d^T  Ki  K2 

-  t4  +  -  =0  [9] 

dx2  k  5  k  5 


Note  that  depends  on  the  thermal  property  (emissivity) 
whereas  K2  depends  on  both  a  thermal  property  (absorptivity) 
and  an  external  factor  (external  heat  source) .  As 
previously  discussed,  external  heat  includes  not  only 
radiation  incident  from  outside  heat  sources,  but  also 
radiation  from  the  fin  that  is  intercepted  by  nearby 
structures.  In  this  sense,  K2  can  be  considered  an 
environmental  parameter.  For  the  free-space  case,  K2  =  0 . 

From  Figure  9  and  the  assumption  (no  heat  transfer 
through  the  tip  of  the  fin) ,  the  boundary  conditions  are 
seen  to  be 

T(0)  =  To 


dT 

dx 


=  0 


x=b 


[10] 
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C.  TRAPEZOIDAL  FIN 
1 .  Geometry 

The  geometry,  terminology  and  coordinate  system  for 
the  longitudinal  fin  of  trapezoidal  profile  are  shown  in 
Figure  10. 


Here: 


Q  =  amount  of  heat  transfer  (W) 

X  =  height  coordinate  along  the  fin  with  the  positive 
orientation  from  base  to  tip  (m) 

dx  =  differential  element  (m) 

T(x)  =  fin  temperature  at  x  (°K) 

Tq  =  temperature  at  fin  base  (°K) 

Tg  =  temperature  at  fin  tip  {°K) 

Ts  =  temperature  of  environment  (for  free  space, 

Ts  =  O^K) 
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L  =  length  of  fin  (in) 
b  =  height  of  fin  (m) 

6(x)  =  thickness  or  width  of  fin  at  x  (m) 

6o  =  thickness  or  width  of  fin  at  the  base  (m) 

6£  =  thickness  or  width  of  fin  at  the  tip  (m) 

2 .  Assumptions 

The  assumptions  are  the  seune  as  those  for  the 
longitudinal  fin  of  rectangular  profile. 

3.  Analysis 

The  analysis  here  is  similar  to  that  for  the 
longitudinal  fin  of  rectangular  profile.  It  differs  in  that 
the  fin  thickness  changes  as  a  function  of  x; 

(6o  -  6e) 

5(x)  =  6g  +  -  X  [11] 


For  each  differential  element  of  the  fin,  the 
difference  between  the  heat  entering  and  leaving  by 
conduction  is 


dQ 


d 

k  L  — 
dx 


[12] 


the  heat  leaving  by  radiation  is 

dQp  =  2  a  e  L  T^  dx 


[13] 
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and  the  heat  incident  from  external  sources  is 


dQg  =  E  a  L  dx 

Application  of  the  general  energy  balance  [ 
differential  element  in  steady-state  establishes 


dQ  »=  dQp  -  dQg 


Substitution  of  equations  [12]  -  [14]  into  equation 
produces 


d 

k  L  — 
dx 


dx 


2  a  e  L  dx  -  E  a  L  dx 


Then  after  substitution  of  *  2  a  e  and  K2  =  E  a 
equation  [16],  simplification  results  in 


d 

dx 


=  Ki  t4  -  K2 


Differentiation  of  the  left  hand  term  gives 


k 


d^T  d5  dT 

6  -  + - 

dx^  dx  dx 


Ki  t4  -  K2 


and  then  further  simplification  provides 

d^T  1  d6  dT  Ki  K2 

—  + - _  -  t4  +  -  -  0 

dx^  6  dx  dx  k  6  k  6 


[14] 
]  to  the 

[15] 

[15] 

[16] 

into 

[17] 

[18] 

[19] 
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From  Figure  10  and  the  assumptions,  the  boundary 


conditions  are  seen  to  be 


T(0)  =  To 


dT 

dx 


x=b 


0 


[20] 


D.  TRIANGULAR  FIN 

The  tip  thickness  for  the  longitudinal  fin  of 
trapezoidal  profile  is  seen  to  be: 

(5o  -  65)  (60  -  6e) 

5(x)  =  5g  +  X  =  6g  +  b  =  6g  [21] 

b  b 

The  longitudinal  fin  of  triangular  profile  is  the 
limiting  case  of  the  trapezoidal  profile.  The  triangular 
fin  has  a  theoretically  "razor  sharp"  tip  with  zero 
thickness  (6g  =  0).  This  results  in  a  singularity  at  6g  =  0 
in  Equation  19.  There  are  several  methods  available  for 
handling  such  a  singularity  [Ref.  11;  p.  48]; 

1.  Integration  by  parts. 

2.  Removal  of  the  singularity  by  addition/ subtraction. 

3.  Integration  through  a  uniform  approach  to  the  limit. 

However,  because  machining  processes  make  it  impossible 

to  achieve  a  fin  tip  with  zero  thickness,  the  elimination  of 
the  singularity  by  one  of  the  foregoing  methods  will  not  be 
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attempted.  The  approximation 


6e  =  0.01  X  6o  [22] 

is  more  realistic  from  the  standpoint  of  fin  construction 
and  also  serves  to  eliminate  the  singularity  at  x  *  0. 


21 


III.  SOLUTION 


A.  PRELIMINARY 
1 .  Review 

The  differential  equations  under  consideration  with 
their  boundary  conditions  are: 

For  the  rectangular  profile: 


d2T 

.  K2 

-  -  - 

T«  +  - 

-  =0 

[23] 

dx2  k  5 

k  6 

dT 

T(0)  =  To 

} 

=  0 

[24] 

dx 

x=b 

For  the  trapezoidal  or  triangular  profiles: 

d2T  1  d6  dT 

Kl 

'  '4-  “  ■  ~  — 

1 1 
II 

o 

[25] 

dx2  6  dx  dx 

k  6 

k  6 

dT 

T(0)  =  To 

7 

=  0 

[26] 

dx 

x=b 

Both  Equations  [23]  and  [25]  are  ordinary,  second 
order,  nonlinear  differential  equations.  Equation  [25], 
however,  possesses  variable  coefficients.  The  presence  of 
the  T^  term  leads  to  the  nonlinearity.  The  requirement  to 
satisfy  conditions  at  two  different  boundaries  classifies 
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Equations  [23]  and  [25]  as  boundary  value  problems.  For  a 
unique  solution,  a  second  order  differential  equation 
requires  two  boundary  conditions.  These  types  of  problems 
are  not  easily  solved  by  normal  analytical  methods  but  can 
be  subjected  to  numerical  procedures  for  solution. 

2.  Second  Order  Differential  Equation 

One  procedure  for  solving  a  second  order  ordinary 
differential  equation  is  to  reduce  it  to  a  system  of  two 
first  order  differential  equations  (canonical  form)  by  an 
appropriate  substitution  and  then  apply  the  techniques 
available  for  solving  first  order  systems  [Ref.  13:  p.  379]. 
By  letting  Y  =  dT/dx  in  Equations  [23]  and  [25],  the  second 
order  ordinary  differential  equations  are  reduced  to  two 
first  order  systems  : 

dT 

—  =  Y 
dx 

[27] 

dY  Ki  K2 

—  =  -  T* - 

dx  k  5  k  6 

with  boundary  conditions: 


T(0)  =  To 
Y(b)  =  0 


[28] 
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and 


dT 

—  =  Y 
dx 

[29] 

dY  1  d6  dT  Ki  ^  K2 
dx  6  dx  dx  k  6  k  6 


with  boundary  conditions: 


T(0)  =  To 
Y(0)  =  0 


[30] 


3.  Boundary  Value  Problem 

An  initial  value  problem  is  a  problem  where  known 
function  and  derivative  values  are  given  at  the  scune  point 
or  boundary  location.  A  boundary  value  problem  is  a 
problem  where  known  function  and  derivative  values  are  given 
at  different  boundary  points  or  locations.  There  are 
numerous  methods  that  can  be  used  to  solve  initial  value 
problems  [Ref.  11:  pp.  50  -  134].  These  include: 

1 .  Taylor  Series  Method 

2.  Euler's  Method 

3 .  Runge-Kutta  Method 

4 .  Multistep  Method 

5.  Predictor-Corrector  Method 

6 .  Adams-Moulton  Method 
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However,  for  boundary  value  problems  there  are 
basically  two  numerical  methods  [Ref.  11:  p.  105]: 

1.  Finite  Difference  Method 

2 .  Shooting  Method 

The  finite  difference  method  uses  a  set  of 
difference  equations  to  approximate  the  solution  to  the 
differential  equation.  The  advantages  to  the  finite 
difference  method  are  its  speed  and  economy  of  memory 
(solving  a  tridiagonal  system  of  equations).  However,  this 
method  runs  into  some  difficulty  when  the  differential 
equation  is  nonlinear  because  it  results  in  a  system  of 
nonlinear  finite  difference  equations  that  must  be  solved  by 
an  iterative  technique.  Thus  the  advantages  of  speed  and 
economy  of  memory  are  largely  eliminated. 

The  shooting  method  guesses  the  slope  of  the 
function  at  one  end  and  then  uses  a  standard  integration 
scheme  to  solve  the  initial  value  problem  to  match  the 
boundary  condition  at  the  other  end  (Figure  11)  [Ref.  13:  p. 
412].  This  procedure  is  repeated  until  the  assumed  solution 
meets  the  specified  or  known  conditions  at  the  boundaries. 
The  advantage  of  the  shooting  method  is  its  simplicity  and 
its  use  of  well  established  methods  to  solve  initial  value 
problems.  Its  disadvantage  is  that  it  is  computationally 
intensive  (dependent  on  the  integration  method  used)  and  may 
take  many  guesses  before  meeting  tolerance. 
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Figure  11  Shooting  Method 


B.  THE  ANALYSIS  PROBLEM 

For  the  analysis  problem,  the  inputs  are  the  fin 
parameters  while  the  output  is  the  eunount  of  heat 
dissipated,  Q.  The  boundary  conditions  are  shown  in  Table 
IV. 

TABLE  IV  ANALYSIS  PROBLEM  BOUNDARY  CONDITIONS 


Left 

Right 

T(0)  =  To 

T(b)  =  ?  =  Te 

dT 

dT 

—  =  ? 

—  =0 

dx  x=0 

dx  x=b 

where;  x  =  0  at  the  base  and  x  =  b  at  the  tip 

Because  the  fin  is  in  steady-state  (there  is  no  heat 
storage),  the  heat  dissipated  must  equal  the  heat  entering 
the  base: 

dT 

Q  =  -  k  A  —  [31] 

dx  x=0 
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Thus,  if  dT/dx  is  known  at  the  base  (x  =  0),  Q  can  be 
calculated.  The  procedure  is  to  guess  dT/dx  at  the  base  (x 
=  0)  and  use  the  shooting  method  to  match  the  other  known 
boundary  condition  dT/dx  =  0  at  the  edge  (x  =  b) .  This 
correct  guess  for  dT/dx  can  be  then  substituted  into 
Equation  [31]  to  yield  the  heat  dissipated,  Q. 

C.  THE  SYNTHESIS  PROBLEM 

For  the  synthesis  problem,  the  input  is  the  amount  of 
heat  to  be  dissipated  while  the  output  is  the  fin  height 
required  for  its  dissipation.  Because  the  fin  is  in 
steady-state  (there  j s  no  heat  storage),  the  heat  dissipated 
must  equal  the  ne^c  entering  the  base: 

dT 

Q  =  -  k  A  —  [32] 

dx  x=0 

However,  unlike  the  analysis  case,  Q  is  known  so  dT/dx 
at  X  =  0  can  be  evaluated  as 

[33] 

where  x  =  0  at  the  base  and  x  =  b  at  the  edge. 

Thus,  in  this  case,  the  boundary  conditions  (see  Figures 
9  and  10)  become 


dT  k  A 

dx  x=0  Q 
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TABLE  V  SYNTHESIS  PROBLEM  BOUNDARY  CONDITIONS 


Left 

Right 

T(0)  =  To 

T{b)  =  ?  =  Tg 

dT 

k  A 

dT 

»  0 

dx 

x=0  Q 

dx 

x=b 

Unfortunately,  b,  the  height  of  the  fin,  which  also 
provides  the  upper  limit  for  integration,  is  unknown.  The 
procedure  here  is  to  shoot  from  the  known  boundary  condition 
dT/dx  =  -kA/Q  at  X  =  0  over  the  guessed  interval  ( fin 
height,  b)  and  match  the  other  known  boundary  condition 
dT/dx  =  0  at  X  =  b. 

D.  ALGORITHM 

1.  Preliminary 

For  the  analysis  problem,  the  procedure  can  be 
viewed  as  the  problem  of  finding  the  root  of  the  function 
w  =  F(2)  where  the  independent  variable,  z,  is  the  guess  of 
the  initial  slope  (dT/dx  at  x  =  0)  and  the  dependent 
variable,  w  =  F(z),  is  the  slope,  dT/dx  at  x  =  b. 

A  root  finding  method  can  then  be  applied  to  find 
the  initial  slope  that  will  match  the  other  boundary 
condition.  Figure  12  illustrates  this  relationship. 

In  similar  fashion,  for  the  synthesis  problem,  the 
procedure  can  be  viewed  as  the  problem  of  finding  the  root 
of  the  function  w  =  F(z)  where  the  independent  variable,  z. 
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Figure  12  Analysis  Problem  as  a  Root  Finding  Problem 

is  the  guess  at  the  height  of  the  fin,  b,  and  the  dependent 
variable,  F(z)  is  once  again  the  final  slope,  dT/dx  at  x  = 
b,  at  the  other  end.  Then  a  root  finding  method  can  again 
be  applied  to  find  the  height  that  will  give  the  correct 
boundary  condition  at  the  fin  tip.  Figure  13  illustrates 
this  relationship. 


Figure  13  Synthesis  Problem  as  a  Root  Finding  Problem 
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2.  Modified  Linear  Interpolation 


There  are  many  root  finding  methods  available  for 
finding  roots  of  nonlinear  equations  [Ref.  14;  p.  27]: 

1 .  Bisection  Method 

2 .  Linear  Interpolation  Method 

3.  Modified  Linear  Interpolation  Method 

4 .  Secant  Method 

5 .  Newton ' s  Method 

6.  Fixed  Point  Iteration  Method 

7.  Muller's  Method 

The  root  finding  method  used  here  is  the  modified 
linear  interpolation  method  (Figure  14).  It  offers  a  good 
combination  of  accuracy  and  speed  in  finding  the  root  of  an 
equation. 


Figure  14  Modified  Linear  Interpolation 


30 


The  algorithm  is  described  in  pseudocode  as  shown  in 
Figure  15  [Ref.  13:  p.  11]. 

Let  X]^  and  X2  bracket  the  root;  ie.,  f(X]^)  and 
f(X2)  are  opposite  in  sign 

Set  SAVE  =  f(xi)}  set  Fl  =  f(X2);  set  F2  =  f(X2) 

Repeat 

X2  “  *1 

Set  x-3  =  XT  -  F2  - 

F2  -  Fl 

If  f(X3)  of  opposite  sign  to  Fl 
Set  X2  =  X3;  set  F2  =  f(X3) 

If  f(X3)  of  same  sign  as  SAVE 
Set  Fl  =  Fl/2 
Endif 

Else 

Set  x^  =  X3;  set  Fl  =  f(X3) 

If  f(X3)  of  same  sign  as  SAVE 
Set  F2  =  F2/2 
Endif 

Endif 

Set  SAVE  =  f(X3) 

Until  |xi  -  X2I  <  XTOL  or  |f(X3)|  <  FTOL 
Figure  15  Modified  Linear  Interpolation  Pseudocode 
where : 

xi  =  X  value  at  left  end  point 
X2  =  X  value  at  right  end  point 
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X3  =  X  value  at  intermediate  point 
f{X2)  =  function  value  at  left  end  point 
f(x2)  =  function  value  at  right  end  point 
f(X3)  =  function  value  at  intermediate  point 
Fl  =  saved  value  of  f(xi) 

F2  =  saved  value  of  f(X2) 

SAVE  =  function  used  for  next  sign  comparison 
XTOL  =  tolerance  in  independent  variable,  x 
FTOL  =  tolerance  in  dependent  variable,  f(x) 

Note  that  halving  the  functional  evaluations  (Fl/2 
or  F2/2)  alleviates  the  difficulty  encountered  by  regular 
linear  interpolation  on  one-sided  approaches  to  the  root  as 
shown  in  Figure  16. 


Figure  16  One-sided  Approach  to  Root 
One  of  the  major  difficulties  in  executing  a  root 
solver  is  to  find  an  interval  that  brackets  the  root 
(bisection,  linear  interpolation,  modified  linear 
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interpolation,  or  secant)  or  that  is  reasonably  close  to  the 
root  (Newton,  fixed  point  iteration,  or  Muller) .  The  method 
used  here  is  a  simple  linear  search  from  the  origin  as  shown 
in  Figure  17  [Ref.  15;  p.  42]. 


f(x) 

1 

interval  1  interval  1  interval  I  .  . 

» 

1  : 
1 

) 

\ 

1 

1 

^  1 

1 

^  1 

_ , _ J 

1 

Figure  17  Search  for  Interval  to  Bracket  Root 
For  the  synthesis  problem,  the  independent  variable 


is  z  =  -dT/dx  and  so  the  algorithm  searches  successive 
intervals  [0,-1],  [-1,-2],  ...  for  functional  values  that 
will  be  opposite  in  sign.  Here,  the  independent  variable 
must  be  negative  to  ensure  that  Q,  the  eunount  of  heat,  is  by 
convention  a  positive  quantity.  The  unit  interval  search  is 
a  compromise  between  too  small  an  interval  that  would 
require  too  many  searches  and  too  large  an  interval  which 


may  miss  the  root  (Figure  18). 
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3 .  Runge-Kutta-Fehlberg 

The  shooting  method  requires  an  integrating  scheme 
to  "march"  or  "shoot"  from  one  boundary  condition  to  the 
other.  As  mentioned  previously,  there  are  many  techniques 
of  numerical  integration: 

1 .  Taylor  Series  Method 

2.  Euler's  Method 

3.  Modified  Euler's  Method 

4 .  Runge-Kutta  Method 
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5. 


Milne  Method 


6  .  Adeuns- Moult  on  Method 

Runoe-Kutta-Fehlberq  offers  the  best  combination  of 
accuracy,  stability  and  adaptability  as  indicated  in  Table 
VI. 


TABLE  VI  COMPARISON  OF  METHODS  FOR  SOLVING  DIFFERENTIAL 

EQUATIONS 


Method 

Global 

Error 

Function 
Evaluations 
Per  Step 

Stability 

Ease  of 
Step  Size 
Adjustment 

Modified  Euler 

0(h2) 

2 

G 

G 

4TH-order  R-K 

0(h5) 

4 

G 

G 

R-K-Fehlberg 

0{h6) 

6 

G 

G 

Milne 

0(h5) 

2 

P 

P 

Adams -Moulton 

0(h5) 

2 

G 

P 

where : 

R-K  =  Runge-Kutta 
h  =  step  size 

0(h2)  =  error  on  the  order  of  h^ 

G  =  good 
P  =  poor 

The  Runge-Kutta-Fehlberg  method  uses  a  4TH-order 
Runge-Kutta  method  to  produce  one  estimate  (Y'n+i)  uses 

a  5TH-order  Runge-Kutta  method  to  produce  another  estimate 
(yN+i)«  It  then  computes  the  error  E  =  y^+i  -  Y'n+I 
adjusts  the  step,  h,  accordingly  while  using  the  y^+i  as  the 
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next  estimate.  Such  a  scheme  is  called  adaptive  or  smart 
because  it  capable  of  estimating  the  error  at  each  step  and 
making  a  corresponding  adjustment  in  step  size  to  remain 
within  the  tolerance.  Because  the  k's  are  the  seune  for  both 
estimates,  only  six  functional  evaluations  are  needed.  The 
power  of  the  method  is  its  high  accuracy,  0(h5),  and 
adaptability.  A  disadvantage  is  that  it  requires  more 
functional  evaluations  than  other  methods.  The  coefficients 
are  difficult  to  derive  algebraically  but  can  be  thought  of 
as  weights  applied  to  previous  estimates  that  are  added  to 
produce  a  future  estimate.  These  values  are  listed  in  Table 
VII. 

TABLE  VII  RUNGE-KUTTA-FEHLBERG  COEFFICIENTS 

ki  =  h  f(XN,yN) 

k2  =  h  f(XN  +  h/4,yN  +  ki/4) 

k3  =  h  f(XN  +  3h/8,yN  +  3ki/32  +  9k2/32) 

k4  =  h  f{XN  +  12h/13,yN  +  1932ki/2197  -  7200k2/2197 
+  7296k3/2197) 

ks  =  h  f(XN  +  h,yN  +  439ki/216  -  8k2  +  3680k3/513 
-  845k4/4104) 

kc  =  h  f(XN  +  h/2,yN  -  8ki/27  +  2k2  -  3544k3/2565 
+  1859k4/4104  -  llk5/40) 

y'N+1  =  YN  +  (25ki/216  +  1408k3/2565  +  2197k4/4104 
-  k5/5)  with  global  error  0(h4) 

YN+l  “  YN  +  (16ki/135  +  6656k3/12825  +  28561k4/54430 
-  9k5/50  +  2kg /55)  with  global  error  0(h5) 

E  =  Yn+1  -  Y'n+1  ®  ki/360  -  128k3/4275  -  2197k4/75240 
+  kg/SO  +  2k5755 


where : 


h  =  step  size 

xjj  =  current  independent  variable,  x 

yjj  =  current  dependent  variable,  y 

f  =  derivative  function  ®  f(x,y) 

-  k5  =  weighted  values  of  f  along  interval 

y'u+i  =  predicted  value  of  y  using  4TH  order  Runge-Kutta 
integration  scheme 

yjj+i  =  predicted  value  of  y  using  5TH  order  Runge-Kutta 
integration  scheme 

0(^4)  =  error  on  the  order  of  h4 

0(h5)  =  error  on  the  order  of  h5 

E  =  error  between  4TH  and  5Th  order  predicted  values 
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IV.  COMPUTER  PROGRAM 


A.  PRELIMINARY 

1.  Structure  Design 

The  computer  progr2un  design  follows  the  concept  of 
structured  design.  Meilir  Page-Jones  in  his  handbook  on 
structure  design  describes  it  as  "the  development  of  a 
blueprint  of  a  computer  system  solution  to  a  problem  that 
has  the  same  components  and  interrelationships  among  the 
components  as  the  original  problem"  [Ref.  16:  p.  3]. 

The  structure  design  concept  has  seven  major  steps 
[Ref.  16:  p.  20]: 

1 .  Problem  recognition 

2 .  Feasibility  Study 

3.  Analysis 

4 .  Design 

5 .  Implementation 

6 .  Testing 

7 .  Maintenance 

The  previous  sections  dealt  with  steps  one  and 
three.  Steps  two,  five  and  seven  refer  to  more  intensive 
software  design  projects  and  are  not  applicable  here.  The 
remaining  sections  will  concentrate  on  steps  four  and  six. 
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2 .  Problem  Review 


Recall  that  the  original  problem  consisted  of  the 
following  elements; 

1.  Environments:  free  space,  non-free  space 

2.  Fin  profiles:  rectangular,  trapezoid,  triangular 

3.  Problem  types:  analysis,  synthesis 

Non-free  space  includes  external  heat  sources.  The 
analysis  problem  determines  the  amount  of  heat  dissipated  by 
a  fin  of  given  dimensions.  The  synthesis  problem  determines 
the  size  of  the  fin  to  dissipate  a  given  quantity  of  heat. 

3 .  Problem  Reduction 

Because  the  free  space  case  is  just  a  special  case 
of  the  non-free  space,  its  solution  will  be  included  in  the 
more  general  case  (K2  =  E  a  =  0  in  Equations  [23]  and  [25]). 

With  this  simplification,  the  problem  can  be 
subdivided  into  six  "sub-problems"  based  on  geometry, 
differential  equation  and  input /outputs  as  listed  in  TABLE 
VIII.  Although  the  triangular  fin  is  a  special  case  of  the 
trapezoidal  fin  (with  tip  thickness  equal  to  1/lOOth  of  base 
thickness),  it  is  still  handled  separately  to  take  advantage 
of  progrcim  modularity  which  is  discussed  later. 

4.  Program  Inputs/Outputs 

In  view  of  the  input /output  format  of  computer 
subroutines,  the  sub-problems  can  also  be  analyzed  from  the 
standpoint  of  inputs/outputs  as  listed  in  TABLE  IX. 
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TABLE  VIII  SUB-PROBLEM  TYPES 


Nr 

Profile 

Problem 

Rectangle 

Analysis 

Rectangle 

Synthesis 

Trapezoid 

Analysis 

Trapezoid 

Synthesis 

5 

Triangle 

Analysis 

6 

Triangle 

Synthesis 

TABLE  IX  INPUT/OUTPUT  RELATIONSHIPS 


Nr 

Shape 

Type 

Input 

Output 

Rectangle 

Anal 

Rectangle 

Syn 

^  r ^ 

Trapezoid 

Anal 

Trapezoid 

Syn 

n 

5 

Triangle 

Anal 

Q » "^E  f  ^ 

6 

Triangle 

Syn 

TQ,L,6Q,k,a,e,E,Q 

^  f  ^E 1 ^ 

where; 

Anal  -  analysis 

Syn  -  synthesis 

Tq  -•  temperature  at  fin  base 

L  -  fin  length 

6  -  fin  width  or  thickness  (rectangle) 

6o  -  fin  width  or  thickness  at  base  (trapezoid  or 
triangle ) 

dg  -  fin  width  or  thickness  at  tip  (trapezoid) 
k  -  conductivity  of  fin 
a  -  absorptivity  of  fin 
6  -  emissivity  of  fin 
Q  -  real  heat  dissipated  by  the  fin 
b  -  fin  height 
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E  -  external  heat  flux 
Tg  -  temperature  at  fin  tip 
n  -  fin  efficiency 

B .  FLOWCHART 

The  flow  chart  logic  is  shown  in  Figure  19.  The  first 
step  is  to  perform  housekeeping  chores  such  as  data 
initialization  and  file  access.  Next  the  units,  fin  profile 
and  problem  type  are  selected  which  in  turn  specifies  the 
fin  parameters  to  be  input.  The  program  then  branches  to 
the  subroutine  that  solves  one  of  the  above  six 
sub-problems.  Each  subprogram  finds  the  root  of  equation 
involving  one  of  the  desired  outputs.  This  requires  solving 
the  associated  second  order,  nonlinear  differential  equation 
describing  the  heat  transfer  under  the  specified  boundary 
conditions.  The  output  is  then  summarized  and  presented. 

The  same  problem  can  be  repeated  with  different  outputs  or  a 
new  problem  can  be  entered  or  the  prograun  can  be  terminated. 
The  abbreviations  used  in  Figure  19  are: 

RECT  =  rectangle 
TRAP  =  trapezoid 
TRI  =  triangle 
ANAL  =  analysis 
SYN  =  synthesis 
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C .  MODULES 


The  basic  principals  of  the  design  step  are; 

1.  Partition  the  system  into  modules. 

2.  Organize  the  modules  into  a  hierarchal  structure. 

The  principal  tool  used  to  do  this  is  the  structure 

chart .  The  structure  chart  "illustrates  the  partitioning  of 
a  system  into  modules  -  showing  their  hierarchy, 
organization,  and  communication"  [Ref.  12;  p.  10]. 

Applying  this  concept  to  Tables  VIII  and  IX,  the 
structure  chart  takes  the  form  as  shown  in  Figure  20.  The 
chart  shows  the  basic  modules  and  their  relationships.  The 
program  for  fin  analysis  is  divided  into  multi-level  modules 
whose  attributes  are  given  in  Table  X.  The  level  refers  to 
the  program  hierarchy  as  shown  in  Table  XI. 

D.  SPECIFICATIONS 

The  program  is  written  in  FORTRAN- 7 7  for  IBM-compatible 
computers.  It  uses  only  very  basic  FORTRAN  statements 
common  to  most  FORTRAN  compilers  to  enhance  compatibility 
with  a  wide  variety  of  computer  environments.  The  program 
has  been  compiled  and  executed  in  the  PC  environment  with 
Microsoft's  Fortran  Compiler.  It  has  also  been  run  on  the 
mainframe  (IBM  360/370)  with  the  WF77  and  FORTVS  compilers. 

The  progreun  is  written  in  double  precision  to  minimize 
round-off  and  truncation  errors,  individually  and  globally, 
that  arise  in  the  computationally  intensive  routines  (root 
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READER  f— IMITIAL  r-INTROOl/CTIOM  ^IffPUT  r-GO  TO  SUBROUTINES  r-OUTPUT  r“PROBLE*l  STATUS  FORMAT 


Figure  20  Computer  Program  Modulea 


TABLE  X  MODULES 


Module 

Level 

Action 

HEADER 

1 

Provides  program  documentation 

INITIAL 

i 

Decla^e^  Variables 

Sets  rn^tial  values 

Opens  fxle  for  data 

Selects  monitor  modes 

INTRO 

1 

Presents  info  menus: 

-  overview 

-  envi^'onments 

-  profiles 

-  problem  types 

-  units 

INPUT 

1 

Processes  input  for: 

-  units 

-  profile 

-  problem  type 

-  tin  parameters 

Summarizes  data  input 

Corrects  data  input 

GOTO 

1 

Based  on  inputs  chooses 
one  of  six  subroutines 

Performs  conversions  as  required 

RTAN 

2 

Solves  rectangle,  analysis  problem 

RTSY 

2 

Solves  rectangle,  synthesis  problem 

TPAN 

2 

Solves  trapezoid,  analysis  problem 

TPSY 

2 

Solves  trapezoid,  synthesis  problem 

TRAN 

2 

Solves  triangle,  analysis  problem 

TRSY 

2 

Solves  triangle,  synthesis  problem 

OUTPUT 

1 

Prints  outputs  to  screen 

Writes  outputs  to  FIN. DAT  file 

PGMSTAT 

1 

Performs  conversion  as  requested 

Sets  up  scune  problem  with 
different  inputs  specified  by  user 
Sets  up  new  problem 

Closes  and  saves  data  file 

Exits  program 

FORMAT 

1 

Input /output  format 

SI 

2 

Converts  to  SI  units 

ENG 

2 

Converts  to  English  units 

FCN 

4 

Computes  function  values 

RKFSY 

5 

Solves  system  of  differential 
equations  by  Runge-Kutta-Fehlberg 

DERIV 

6 

Defines  derivatives 

MDLIN 

3 

Solves  for  root  of  function  by 
modified  linear  interpolation 
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TABLE  XI  PROGRAM  HIERARCHY 


solver  and  integrator) .  Most  inputs  will  have  at  best  4-6 
significant  figures  due  to  physical  measurement  limitations, 
so  the  accuracy  of  the  results  can  not  be  expected  to  exceed 
that . 

The  progrcun  is  completely  independent  and  invokes  no 
external  library  routines.  All  numerical  subroutines  for 
root  finding  and  integration  have  been  incorporated  into  the 
progreim. 

The  simplicity  mentioned  above  has  been  achieved  at  the 
cost  of  program  length  -  nearly  4500  lines  of  code. 
Additional  length  also  resulted  from  the  duplication  of  some 
sections  to  maintain  program  modularity.  For  development 
and  subsequent  improvement,  it  is  important  that  individual 
modules  have  the  capability  of  being  removed,  modified  and 
reinserted  with  minimum  interference.  To  this  end,  all  main 
program  elements  and  sub-elements  can  be  easily  removed  and 
independently  operated  with  simple  "drivers". 

Program  specifications  are  summarized  in  Table  XII. 


46 


TABLE  XII  PROGRAM  SPL3IFICATIONS 


Parameter 

Specification 

Language 

FORTRAN  77 

Lines  of  code 

4425 

Memory :  Fortran  code 

140020 

Object  code 

102455 

Execution  code 

107370 

Total : 

349845 

Compiler 

Microsoft  FORTRAN  V5.0 

Compile  time 

3.5  minutes 

Computer 

IBM  clone  ( 386/387/16 . OMhz ) 

Execution  time 

Variable  (see  examples) 

Precision 

Double 

External  programs/ 

None 

libraries  required 

E.  MENUS 

The  program  is  executed  in  the  interactive  mode.  The 
user  interfaces  with  the  progreim  through  a  series  of  menus. 
The  two  interfacing  modes  are: 

1 .  Line  edit 

2 .  Screen  edit 

In  the  line  edit  mode,  the  user  interacts  with  the 
program's  commands  as  they  are  "scrolled"  to  the  screen. 
This  enables  the  last  few  commands  to  always  be  on  the 
screen.  In  the  screen  edit  mode,  there  is  one  command  per 
screen  and  the  screen  is  cleared  before  the  next  command, 
thereby  minimizing  clutter  on  the  screen.  The  screen  edit 
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mode  can  only  by  used  on  a  typical  24  line  PC  monitor.  It 
will  not  "synchronize*  with  other  monitors  or  with  the  video 
terminals  associated  with  the  mainfr2une.  The  line  edit  mode 
should  be  used  in  these  cases. 

The  menus  are  straightforward  and  consist  of  simple 
questions  or  commands  to  enter  inputs.  Table  XIII  is  a 
listing  of  the  menus. 


TABLE  XIII  MENUS 


NR 

Query/ Input 

1 

"CAPS  LOCK"  selection 

2 

Viewing  mode  selection 

3 

Introduction 

4 

Units  selection 

5 

Fin  profile  selection 

6 

Problem  type  selection 

7 

Base  temperature  input 

8 

Length  input 

9 

Thickness  input 

10 

Conductivity  input 

11 

Absorptivity  input 

12 

Emissivity  input 

13 

External  heat  input 

14 

Heat  entering  base  (synthesis)  input 

15 

Height  (analysis)  input 

16 

Incorrect  entry  input 

17 

Output  conversion 

18 

Problem  selection:  same  problem  with  changes 

19 

Problem  selection:  new  problem 

20 

Quit 

F .  ERROR  ROUTINES 

Several  error  routines  have  been  built  into  the  progrsun. 
These  can  be  roughly  divided  into  three  major  categories: 
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1.  Input  error  detection 

2 .  Output  error  detection 

3.  Progreun  error  detection 

The  first  category  concerns  errors  made  by  the  operator 
during  the  course  of  menu  interrogation.  Although  no 
program  can  be  completely  "idiot  proof",  an  attempt  has  been 
made  to  minimize  input  errors  by  restricting  the  input  range 
to  realistic  values  as  shown  in  Table  XIV.  Any  input 
outside  this  range  will  result  in  an  "INCORRECT  ENTRY"  flag 
along  with  a  repetition  of  the  menu  command. 


TABLE  XIV  INPUT  PARAMETER  RANGE 


Parameter 

Range 

y  or  N 

Capital  "Y"  or  capital  "N" 

1,2, .. . 

Numbered  response  "1","2",... 

Temp  at  fin  base  (Tq) 

a  -273.150k  (-460. O^F) 

Fin  length  (L) 

>  0 

Base  thickness  (5  or  6q) 

>  0 

Tip  thickness  (fig) 

>  0 

Conductivity  (k) 

>  0 

Absorptivity  (a) 

0  <  a  £  1 

Emissivity  (e) 

0  <  e  £  1 

External  heat  flux  (E) 

i  0 

Heat  input  thru  base  (Q) 

>  0 

Fin  height  (b) 

>  0 

Similarly,  any  values  outside  the  ranges  for  the  outputs 
listed  in  Table  XV  will  result  in  a  flag  to  check  the  inputs 
for  correctness. 

Finally,  there  are  internal  error  checks  to  minimize  the 
possibility  of  a  "hung"  progreun.  Here,  the  term  "hung" 
program  means  that  the  algorithm  is  caught  in  an  "infinite" 
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TABLE  XV  OUTPUT  PARAMETER  RANGE 


Parameter 

Range 

Heat  input  thru  base  (Q) 
Fin  height  (b) 

Temp  at  fin  tip  (Tg) 
Efficiency  (fl) 

>  0 
>  0 

s  -273.15°K,  -460. 0°F 

0  s  n  s  1 

loop;  i.e.,  the  problem  oscillates  or  diverges.  No 
numerical  procedure  can  or  should  be  made  so  "robust"  that 
it  converges  regardless  of  the  inputs.  In  this  sense,  the 
best  defense  is  to  ensure  that  the  inputs  are  not  only 
restricted  to  the  ranges  listed  above,  but  that  they  also 
appear  reasonable.  Fins  with  thicknesses  of  a  hundred 
meters  or  with  base  temperatures  of  ten  thousand  degrees  are 
unreasonable . 

The  three  major  numerical  algorithms  (interval  search, 
root  solver  and  differential  equation  solver)  have  several 
methods  of  stopping.  The  interval  search  routine  is  limited 
to  a  fixed  number  of  intervals  (100)  to  be  searched.  If  a 
root  is  not  found  within  the  first  100  intervals,  the  user 
is  asked  if  he  wants  to  continue  the  search  for  another  100 
intervals.  This  procedure  is  repeated  until  an  interval  is 
found  or  the  user's  patience  is  exhausted.  The  user  may 
also  request  that  the  function  values  be  printed  to  the 
screen.  This  will  give  some  idea  of  the  function's 
behavior.  Recall  that  for  a  function  to  have  a  root  at  x, 
f(x)  must  equal  zero  and  hence  any  interval  that  brackets 
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the  root  must  have  function  values  at  the  endpoints  that  are 
opposite  in  sign.  So  by  examining  the  function  values,  the 
user  can  get  a  "feel"  for  whether  or  not  the  function  values 
are  approaching  zero  with  a  resultant  sign  change.  Note 
that  this  is  not  an  absolute  guarantee  since  a  function  can 
increase  then  decrease  or  vice  versa. 

Both  the  root  and  differential  equation  solvers 
terminate  through  one  of  the  following  mechanisms: 

1.  Successive  estimated  values  are  within  a  specified 
tolerance. 

2.  A  fixed  number  of  iterations  is  exceeded. 

The  termination  criteria  are  listed  in  Table  XVI .  As  is 
the  case  for  the  interval  search  algorithm,  the  user  can 
view  the  function  values  to  get  a  rough  idea  of  convergence 
or  divergence. 

TABLE  XVI  ROOT  SOLVER  AND  DIFFERENTIAL  EQUATION  SOLVER 

STOPPING  CRITERIA 


Routine 

Parameter 

Stopping  criteria 

Modified 

X  tolerance 

0.0001 

linear 

interpolation 

F(X)  tolerance 

0.00001 

Iterations 

50 

Runge-Kutta- 

Y  tolerance 

0.0001 

Fehlberg 

Iterations 

100 
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V.  EXAMPLES 


A.  PRELIMINARY 

The  following  ex£unples  demonstrate  the  operation  of  the 
progrcun.  The  results  compared  favorably  with  "manual"  plots 
[Ref  7,  pp.  242  -  272].  The  input/output  summaries  are 
copies  of  ones  that  were  actually  written  to  the  file 
FIN. DAT.  Note  that  the  synthesis  problems  are  intentionally 
"inverses"  of  the  analysis  problems  so  as  to  further  check 
results.  Also  the  execution  times  on  an  IBM  clone 
(386/387/16.0  Mh2)  are  included  in  parenthesis  for 
comparison.  The  times  are  a  function  of  the  machine  as  well 
as  the  problem  itself  (geometry  and  fin  parameters).  The 
triangular  fin/synthesis  problem  is  particularly  slow 
because  of  the  small  step  used  in  integration. 

B.  RECTANGULAR  FIN 

1.  Analysis  (29  seconds) 

A  rectangular  fin  has  a  base  temperature  of  77  ®C,  a 
length  of  4  meters,  a  thickness  of  0.635  centimeters,  a 
height  of  50  centimeters,  a  conductivity  of  152  Watts/meter, 
an  absorptivity  of  zero,  and  an  emissivity  of  0.85.  There 
are  no  external  heat  sources.  What  is  the  eunount  of  heat 
dissipated  by  the  fin?  What  is  the  temperature  at  the  fin 
tip?  What  is  the  fin  efficiency? 
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INPUT/OUTPUT  SUMMARY; 


**********************  input  *************************** 


(1) 

UNITS 

S 

SI 

(2) 

FIN 

=s 

RECTANGLE 

(3) 

PROBLEM 

= 

ANALYSIS 

(4) 

BASE  TEMPERATURE 

s= 

77.0000 

C 

(5) 

FIN  LENGTH 

3S 

4.00000 

M 

(6) 

BASE  THICKNESS 

S= 

.635000E-02 

M 

(7) 

CONDUCTIVITY  CONSTANT 

= 

152.000 

W/M-C 

(8) 

ABSORPTION  COEFFICIENT= 

.000000 

(9) 

EMISSIVITY 

S 

.850000 

(10) 

EXTERNAL  HEAT  FLUX 

= 

.000000 

W 

(11) 

FIN  HEIGHT 

ss 

.500000 

M 

**********************  OUTPUT 

(12) 

HEAT  INPUT  THRU  BASE 

= 

1510.83 

W 

(13) 

EDGE  TEMPERATURE 

-3.32502 

c 

(14) 

FIN  EFFICIENCY 

= 

.521312 

2.  Synthesis  (43  seconds) 

A  rectangular  fin  has  a  base  temperature  of  77 
length  of  4  meters,  a  thickness  of  0.635  centimeters,  a 
conductivity  of  152  Watts/meter,  an  absorptivity  of  zero, 
and  an  emissivity  of  0.85.  There  are  no  external  heat 
sources.  What  is  the  height  required  for  the  fin  to 
dissipate  1510.83  Watts  of  heat  entering  through  the  base 
(see  the  previous  problem)?  What  is  the  temperature  at  the 
fin  tip?  What  is  the  fin  efficiency? 

INPUT/OUTPUT  SUMMARY: 


**********************  input 

(1)  UNITS 

(2)  FIN 

( 3 )  PROBLEM 

(4)  BASE  TEMPERATURE 

(5)  FIN  LENGTH 

(6)  BASE  THICKNESS 


iritifir’kieir'kitirifitie’k’kififirificitifitieicit 

SI 

RECTANGLE 
SYNTHESIS 
77.0000  C 

4.00000  M 

.635000E-02  M 


(7)  CONDUCTIVITY  CONSTANT  = 

(8)  ABSORPTION  COEFFICIENT= 

(9)  EMISSIVITY 

(10)  EXTERNAL  HEAT  FLUX 

(11)  HEAT  INPUT  THRU  Bi.SE  = 


152.000  W/M-C 

.000000 
.850000 
.000000  W 

1510.83  W 


*****************11****  OUTPUT  ************************** 


(12)  FIN  HEIGHT  =  .500015  M 

(13)  EDGE  TEMPERATURE  =  -3.32527  C 

(14)  FIN  EFFICIENCY  =  .521297 

C.  TRAPEZOIDAL  FIN 

1.  Analysis  (9  seconds) 

A  trapezoidal  fin  has  a  base  temperature  of  167  ®C, 
a  length  of  1.143  meters,  a  base  thickness  of  0.9525 
centimeters,  an  edge  thickness  of  0.47625  centimeters,  a 
height  of  15.24  centimeters,  a  conductivity  of  33 
Watts /meter,  an  absorptivity  of  zero,  and  an  emissivity  of 
0.95,  There  are  no  external  heat  sources.  What  is  the 
amount  of  heat  dissipated  by  the  fin?  What  is  the 
temperature  at  the  fin  tip?  What  is  the  fin  efficiency? 

INPUT/OUTPUT  SUMMARY: 


(1) 

UNITS 

= 

SI 

(2) 

FIN 

= 

TRAPEZOID 

(3) 

PROBLEM 

s 

ANALYSIS 

(4) 

BASE  TEMPERATURE 

s 

167.000 

C 

(5) 

FIN  LENGTH 

— 

1.14300 

M 

(6) 

BASE /EDGE  THICKNESSES 

= 

.952500E-02 

.476250E 

(7) 

CONDUCTIVITY  CONSTANT 

s 

33.0000 

W/M-C 

(8) 

ABSORPTION  COEFFICIENT= 

.000000 

(9) 

EMISSIVITY 

.950000 

(10) 

EXTERNAL  HEAT  FLUX 

= 

.000000 

W 

(11) 

FIN  HEIGHT 

= 

.152400 

M 
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**********************  OUTPUT  ************************** 

(12)  HEAT  INPUT  THRU  BASE  =  412.048  W 

(13)  EDGE  TEMPERATURE  =  77.7946  C 

(14)  FIN  EFFICIENCY  =  .584976 

2.  Synthesis  (11  seconds) 

A  trapezoidal  fin  has  a  base  temperature  of  167  ®C, 
a  length  of  1.143  meters,  a  base  thickness  of  0.9525 
centimeters,  an  edge  thickness  of  0.47625  centimeters,  a 
conductivity  of  33  Watts /meter,  an  absorptivity  of  zero,  and 
an  emissivity  of  0.95.  There  are  no  external  heat  sources. 
What  is  the  height  required  for  the  fin  to  dissipate  412.048 
Watts  of  heat  entering  through  the  base  (see  the  previous 
problem) ?  What  is  the  temperature  at  the  fin  tip?  What  is 
the  fin  efficiency? 

INPUT/OUTPUT  SUMMARY; 

**********************  input  *************************** 


(1) 

UNITS 

= 

SI 

(2) 

FIN 

s 

TRAPEZOID 

(3) 

PROBLEM 

= 

SYNTHESIS 

(4) 

BASE  TEMPERATURE 

= 

167.000 

C 

(5) 

FIN  LENGTH 

1.14300 

M 

(6) 

BASE/EDGE  THICKNESSES 

= 

.952500E-02 

476250E-02 

(7) 

CONDUCTIVITY  CONSTANT 

= 

33.0000 

W/M-C 

(8) 

ABSORPTION  COEFFICIENT= 

.000000 

(9) 

EMISSIVITY 

= 

.950000 

(10) 

EXTERNAL  HEAT  FLUX 

= 

.000000 

W 

(11) 

HEAT  INPUT  THRU  BASE 

412.048 

W 

★  ★  ★  * 

******************  OUTPUT 

Icir'k'k'k'k'kiricicleit'k'k'kir'kitifk’kifirieicic 

(12) 

FIN  HEIGHT 

S 

.152399 

M 

(13) 

EDGE  TEMPERATURE 

s: 

77.7945 

C 

(14) 

FIN  EFFICIENCY 

= 

.584980 
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D.  TRIANGULAR  FIN 

1.  Analysis  (14  seconds) 

A  triangular  fin  has  a  base  temperature  of  167  a 
length  of  1.143  meters,  a  base  thickness  of  0.9525 
centimeters,  a  height  of  16.5  centimeters,  a  conductivity  of 
33  Watts/meter,  an  absorptivity  of  zero,  and  an  emissivity 
of  0.95.  There  are  no  external  heat  sources.  What  is  the 
aunount  of  heat  dissipated  by  the  fin?  What  is  the 
temperature  at  the  fin  tip?  What  is  the  fin  efficiency? 

INPUT /OUTPUT  SUMMARY: 

**********************  input  *************************** 


(1)  UNITS 

(2)  FIN 

( 3 )  PROBLEM 

(4)  BASE  TEMPERATURE 

(5)  FIN  LENGTH 

(6)  BASE  THICKNESS 

(7)  CONDUCTIVITY  CONSTANT  = 

(8)  ABSORPTION  COEFFICIENT* 

(9)  EMISSIVITY 

(10)  EXTERNAL  HEAT  FLUX 

(11)  FIN  HEIGHT 

**********************  OUTPUT 


SI 

TRIANGLE 
ANALYSIS 
167.000  C 

1.14300  M 

.952500E-02  M 
33.0000  W/M-C 

.000000 
.950000 
.000000  W 

.165000  M 

************************** 


(12)  HEAT  INPUT  THRU  BASE  =  400.040  W 

(13)  EDGE  TEMPERATURE  =  40.4865  C 

(14)  FIN  EFFICIENCY  =  .524560 

2.  Synthesis  (3.7  minutes) 

A  triangular  fin  has  a  base  temperature  of  167  ®C,  a 
length  of  1.143  meters,  a  base  thickness  of  0.9525 
centimeters,  a  conductivity  of  33  Watts/meter,  an 
absorptivity  of  zero,  and  an  emissivity  of  0.95.  There  are 
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no  external  heat  sources.  What  is  the  height  required  for 
the  fin  to  dissipate  400.040  Watts  of  heat  entering  through 
the  base  (see  the  previous  problem)?  What  is  the 
temperature  at  the  fin  tip?  What  is  the  fin  efficiency? 

INPUT/OUTPUT  SUMMARY; 

**********************  input  *************************** 


(1)  UNITS 

(2)  FIN 

( 3 )  PROBLEM 

(4)  BASE  TEMPERATURE 

(5)  FIN  LENGTH 

(6)  BASE  THICKNESS 

(7)  CONDUCTIVITY  CONSTANT  = 

(8)  ABSORPTION  COEFFICIENT= 

(9)  EMISSIVITY 

(10)  EXTERNAL  HEAT  FLUX 

(11)  HEAT  INPUT  THRU  BASE  = 


SI 

TRIANGLE 
SYNTHESIS 
167.000  C 

1.14300  M 

.952500E-02  M 
33.0000  W/M-C 

.000000 
.950000 
.000000  W 

400.040  W 


**********************  OUTPUT  ************************** 


(12)  FIN  HEIGHT  =  .165480  M 

(13)  EDGE  TEMPERATURE  =  42.2041  C 

(14)  FIN  EFFICIENCY  =  .523037 


APPENDIX  A 


SI  UNITS 


Neiine 

Ouantitv 

Symbol 

meter 

length 

m 

square  meter 

area 

m^ 

cubic  meter 

volume 

m^ 

kilogram 

mass 

kg 

second 

time 

s 

Newton 

force 

N 

degree  Kelvin 

temperature 

OK 

degree  Celsius 

temperature 

OC 

meter/second 

velocity 

m/s 

newton/square  meter 

pressure 

N/m2 

joule 

work 

J 

joule/second  (watt) 

power 

J/s  (W) 
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APPENDIX  B 


ENGLISH  UNITS 


Ncune 

foot 

square  foot 

cubic  foot 

pound 

second 

pound  force 

degree  Fahrenheit 

degree  Rankine  (absolute) 

foot/second 

ft/square  foot 

British  Thermal  Unit 
foot-pound  force 

British  Thermal  Unit/sec 
horsepower 

foot-pound  force/second 


Svmbol 

length 

ft 

area 

ft2 

volume 

ft3 

mass 

lb 

time 

s 

force 

Ibf 

temperature 

Op 

temperature 

Or 

velocity 

ft/ S 

pressure 

Ibf /m2 

work 

BTU 

work 

ft-lbf 

power 

BTU/s 

power 

hp 

power 

ft-lbf/s 
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APPENDIX  C 


PROGRAM  LISTING 

C*  ***********************  ************  ***■»*■»•******  **************  ******* 

c 

C  INDEX;  USE  THE  THREE  ALPHANUMERIC  CHARACTERS  BELOW  AND  THE  SEARCH/ 

C  FIND  FUNCTION  OF  YOUR  EDITOR  TO  QUICKLY  LOCATE  THE  DESIRED 

C  SECTION. 

C 


c 

CIA 

- 

HEADER 

c 

CIB 

- 

MAIN 

c 

ClC 

- 

INTRC 

>  MENUS 

c 

ClD 

- 

INPUl 

DATA 

c 

ClE 

- 

SUMMARIZE  DATA 

INPUT 

c 

CIF 

- 

GOTO 

SUBROUTINES 

c 

ClG 

- 

SUMMARIZE  DATA 

OUTPUT 

c 

ClH 

- 

FORMAT  STATEMENTS 

c 

C2A 

- 

RTAN 

HEADER 1 

c 

C2B 

- 

RTAN 

MAINl 

c 

C2C 

- 

RTAN 

FCNl 

c 

C2D 

- 

RTAN 

RKFSYl 

c 

C2E 

- 

RTAN 

DERIVl 

c 

C2F 

- 

RTAN 

MDLINl 

c 

C3A 

- 

RTSY 

HEADER2 

c 

C3B 

- 

RTSY 

MAIN2 

c 

C3C 

- 

RTSY 

FCN2 

c 

C3D 

- 

RTSY 

RKFSY2 

c 

C3E 

- 

RTSY 

DERIV2 

c 

C3F 

- 

RTSY 

MDLIN2 

c 

C4A 

- 

TPAN 

HEADER3 

c 

C4B 

- 

TPAN 

MAIN3 

c 

C4C 

- 

TPAN 

FCN3 

c 

C4D 

- 

TPAN 

RKFSY3 

c 

C4E 

- 

TPAN 

DERIV3 

c 

C4F 

- 

TPAN 

MDLIN3 

c 

C5A 

- 

TPSY 

HEADER4 

c 

C5B 

- 

TPSY 

MAIN4 

c 

CSC 

- 

TPSY 

FCN4 

c 

C5D 

- 

TPSY 

RKFSY4 

c 

C5E 

- 

TPSY 

DERIV4 

c 

CSF 

- 

TPSY 

MDLIN4 

c 

C6A 

- 

TRAN 

HEADERS 

c 

C6B 

- 

TRAN 

MAINS 

c 

C6C 

- 

TRAN 

FCNS 

c 

C6D 

- 

TRAN 

RKFSYS 

c 

C6E 

- 

TRAN 

DERIVS 

c 

C6F 

- 

TRAN 

MDLINS 

c 

C7A 

- 

TRSY 

HEADERS 

c 

C7B 

- 

TRSY 

MAIN6 
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c 

C7C  -  TRSY: 

FCN6 

c 

C7D  -  TRSY: 

RKFSY6 

c 

C7E  -  TRSY: 

DERIV6 

c 

C7F  -  TRSY: 

MDLIN6 

c 

C8A  -  SI 

c 

C9A  -  ENG 

ClA  *****  HEADER  ***************************************************** 
C 

C  MAIN  PRCXSRAM:  FIN. FOR 


PURPOSE:  THIS  PROGRAM  SOLVES  THE  FOLLOWING  LONGITUDINAL  FIN 

HEAT  TRANSFER  PROBLEMS: 

FIN  PROFILES:  RECTANGULAR, TRAPEZOIDAL,  &  TRIANGULAR 
PROBLEM  TYPES:  ANALYSIS  &  SYNTHESIS 
ENVIRONMENTS:  FREE  SPACE  &  NON-FREE  SPACE 


C  PROGRAMMER:  D.  R.  JOHNSON 
C 

C  DATE:  10  JUN  90 
C 

C - 

C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c- 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


ANALYSIS  PROBLEM: 

INPUTS : 

RECTANGLE 

BASE  TEMPERATURE 
LENGTH 

BASE  THICKNESS 

CONDUCTIVITY 

ABSORPTIVITY 

EMISSIVITY 

HEIGHT 

EXTERNAL  HEAT  FLUX 


TRAPEZOID 


TRIANGLE 


BASE  TEMPERATURE 
LENGTH 

BASE  THICKNESS 

TIP  THICKNESS 

CONDUCTIVTY 

ABSORPTIVITY 

EMISSIVITY 

HEIGHT 

EXTERNAL  HEAT  FLUX 


BASE  TEMPERATURE 
LENGTH 

BASE  THICKNESS 
* 

CONDUCTIVITY 

ABSORPTIVITY 

EMISSIVITY 

HEIGHT 

EXTERNAL  HEAT  FLUX 


*  TIP/BASE  THICKNESS  RATIO  ASSUMED  =  0.01 
OUTPUTS : 

RECTANGLE  TRAPEZOID  TRIANGLE 


HEAT  DISSIPATED 
TIP  TEMPERATURE 
EFFICIENCY 


HEAT  DISSIPATED 
TIP  TEMPERATURE 
EFFICIENCY 


HEAT  DISSIPATED 
TIP  TEMPERATURE 
EFFICIENCY 
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n  o 


c 


C  SYNTHESIS  PROBLEM: 

C 

C  INPUTS ; 

C 

C  RECTANGLE  TRAPEZOID  TRIANGLE 

C  -  -  - 

C  BASE  TEMPERATURE  BASE  TEMPERATURE  BASE  TEMPERATURE 

C  LENGTH  LENGTH  LENGTH 

C  BASE  THICKNESS  BASE  THICKNESS  BASE  THICKNESS 

C  TIP  THICKNESS  * 

C  CONDUCTIVITY  CONDUCTIVTY  CONDUCTIVITY 

C  ABSORPTIVITY  ABSORPTIVITY  ABSORPTIVITY 

C  EMISSIVITY  EMISSIVITY  EMISSIVITY 

C  EXTERNAL  HEAT  FLUX  EXTERNAL  HEAT  FLUX  EXTERNAL  HEAT  FLUX 

C  HEAT  INPUT  HEAT  INPUT  HEAT  INPUT 

C 

C  *  TIP/BASE  THICKNESS  RATIO  ASSUMED  =  0.01 
C 

C  OUTPUTS : 

C 

C  RECTANGLE  TRAPEZOID  TRIANGLE 

C  -  -  - 

C  HEIGHT  HEIGHT  HEIGHT 

C  TIP  TEMPERATURE  TIP  TEMPERATURE  TIP  TEMPERATURE 

C  EFFICIENCY  EFFICIENCY  EFFICIENCY 

C 

C - 

c 

C  PARAMETERS : 

C 

C  TB  -  TEMPERATURE  AT  THE  BASE  OF  THE  FIN 

C  TE  -  TEMPERATURE  AT  THE  TIP  OF  THE  FIN 

C  HT  -  HEIGHT  OF  THE  FIN 
CL-  LENGTH  OF  THE  FIN 

C  DEL  -  THICKNESS  OF  FIN  ( RECTANGLULAR) 

C  DELO  -  THICKNESS  OF  FIN  AT  BASE  (TRAPEZOIDAL  OR  TRIANGULAR) 

C  DELE  -  THICKNESS  OF  FIN  AT  TIP  (TRAPEZOIDAL) 

C  K  -  CONDUCTIVITY  OF  THE  FIN 
C  K1  -  CONSTANT  =  2  *  SB  *  EMIS 

C  K2  -  CONSTANT  =  E  *  ABS 

C  ABS  -  ABSORPTIVITY  OF  THE  FIN 
C  EMIS  -  EMISSIVITY  OF  THE  FIN 
C  E  -  EXTERNAL  HEAT  INCIDENT  ON  THE  FIN 
C  EFF  -  EFFICIENCY  OF  THE  FIN 
C  Q  -  REAL  HEAT  DISSIPATED  BY  THE  FIN 
C  QI  -  IDEAL  HEAT  DISSIPATED  BY  THE  FIN 
C  SB  -  STEFAN-BOLTZMAN  CONSTANT 
C 

C - 

C 
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C  SUBROUTINES ; 

C 

C  RTAN  -  SOLVES  RECTANGULAR  FIN,  ANALYSIS  PROBLEM 

C  RTSY  -  SOLVES  RECTANGULAR  FIN,  SYNTHESIS  PROBLEM 

C  TPAN  -  SOLVES  TRAPEZOIDAL  FIN,  ANALYSIS  PROBLEM 

C  TPSY  -  SOLVES  TRAPEZOIDAL  FIN,  SYNTHESIS  PROBLEM 

C  TRAN  -  SOLVES  TRIANGULAR  FIN,  ANALYSIS  PROBLEM 

C  TRSY  -  SOLVES  TRIANGULAR  FIN,  SYNTHESIS  PROBLEM 

C 

C - 

C 

C  UNITS; 

C 

C  TYPE  SI  (INTERNATIONAL)  ENGLISH 

C  LINEAR  METERS  FEET 

C  TEMPERATURE  CENTIGRADE  FAHRENHEIT 

C  HEAT  WATTS  BTU/HR 

C 

ClB  *****  MAIN  ******************************************************* 

C 

REAL* 8  TB , TE , HT , L , DEL , DELO , DELE , K, ABS , EMIS , E , Q, QI , EFF 
INTEGER  PASS, FLAG 

CHARACTER* 2  ANS , UNITS , FIN , TYPE, MODE 

DATA  TB, TE, HT, L, DEL, DELO, DELE, K, ABS, EMIS, E,Q,QI, EFF/ 14*0. ODD/ 
FLAG  =  0 


C 

c -  OPEN  FILE  FIN. DAT  TO  HOLD  DATA  FOR  PRINTER  - 

C 

OPEN(UNIT=9, 

&  FILE='FIN.DAT' , 

&  ACCESS= • SEQUENTIAL • , 

&  FORM= • FORMATTED • , 

&  STATUS= • UNKNOWN ’ ) 

C 

ClC  *****  INTRODUCTION  MENUS  ***************************************** 
C 

c -  ENSURE  IN  "CAPS  LOCK"  MODE  - 

C 


1  PRINT  800 
READ  9000, ANS 

IF(ANS.NE. ’Y' .AND.ANS.NE. ‘N’ )  THEN 
PRINT  1200 
GOTO  1 
END  IF 

IF(ANS.EQ. 'N' )  GOTO  1 
C 
C 

C -  SELECT  MONITOR  MODE  - 

C 

2  PRINT  900 
READ  9 000, MODE 

IF(MODE.NE. '1' .AND. MODE. NE. '2' )  THEN 
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PRINT  1200 
GOTO  2 
ENDIF 
C 

C - SKIP  MENUS - 

C 

IF(MODE.EQ. • 1* )  PRINT  1000 
IF(MODE.EQ. •2* )  PRINT  1001 
3  PRINT  1100 
READ  9000, ANS 

IF(AMS.NE. ’ Y' .AND.ANS.NE. 'N* )  THEN 
PRINT  1200 
GOTO  3 
ENDIF 

IF(ANS.EQ. 'Y' )  GO  TO  44 
C 

C - INTRO  MENU  -  GENERAL - 

C 

IF(MODE.EQ. • 1* )  PRINT  1000 
IF(MODE.EQ. *2 • )  PRINT  1001 
PRINT  1300 
PRINT  1301 
10  PRINT  1400 
READ  9000, ANS 
IF(ANS.NE.*  •)  THEN 
PRINT  1200 
GOTO  10 
ENDIF 
C 

C - INTRO  MENU  -  FIN  ANALYSIS  INPUT 

C 

IF(MODE.EC!.  •  1’ )  PRINT  1000 
IF{MODE.EQ. *2 • )  PRINT  1001 
PRINT  1500 
15  PRINT  1400 
READ  9000, ANS 
IF(ANS.NE.'  •)  THEN 
PRINT  1200 
GOTO  15 
ENDIF 
C 

C - INTRO  MENU  -  FIN  ANALYSIS  OUTPUT 

C 

IF(MODE.EQ. • 1' )  PRINT  1000 
IF(MODE.EQ. ’2 • )  PRINT  1001 
PRINT  1600 
20  PRINT  1400 
READ  9000, ANS 
IF(ANS.NE.*  •)  THEN 
PRINT  1200 
GOTO  20 
ENDIF 
C 
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c - INTRO  MENU  -  FIN  SYNTHESIS  INPUT - 

C 

IF(MODE.EQ. • !• )  PRINT  1000 
IF(MODE.EQ. •2’ )  PRINT  1001 
PRINT  1700 
25  PRINT  1400 
READ  9000, ANS 
IF(ANS.NE.*  •)  THEN 
PRINT  1200 
GOTO  25 
ENDIF 
C 

C -  INTRO  MENU  -  FIN  SYNTHESIS  OUTPUT  - 

C 

IF(MODE.EQ. • 1 • )  PRINT  1000 
IF(MODE.EQ. '2 • )  PRINT  1001 
PRINT  1800 
30  PRINT  1400 
READ  9000, ANS 
IF(ANS.NE.’  •)  THEN 
PRINT  1200 
GOTO  30 
ENDIF 
C 

c - INTRO  MENU  -  UNITS  MODE - 

C 

35  IF(MODE.EQ. ' 1 • )  PRINT  1000 
IF(MODE.EQ. '2 ■ )  PRINT  1001 
PRINT  1900 
PRINT  1400 
RE.  D  9000, ANS 
IF(ANS.NE.*  •)  THEN 
PRINT  1200 
GOTO  35 
ENL'IF 
C 

c - INTRO  MENU  -  NOTES - 

C 

IF  MODE.EQ.'l’)  PRINT  1000 
IF(MODE.EQ. '2 ’ )  PRINT  1001 
PP.NT  2000 
PP:NT  2001 
40  PP:NT  1400 
RE,J>  9000, ANS 
IF  ANS.NE. ■  • )  THEN 

''RINT  1200 
GOTO  40 
ENDIF 
C 

ClD  *****  INPUT  DATA  ************************************************* 

C 

C -  SELECT  UNITS  MODE  - 

C 
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44  PASS  =  0 

45  IF(MODE.EQ. • 1* )  PRINT  1000 
IF(MODE.EQ. *2' )  PRINT  1001 

46  PRINT  2100 
READ  9000, UNITS 

IF(UNITS.NE. • 1* .AND. UNITS. NE. ’2* )  THEN 
PRINT  1200 
GOTO  46 
ENDIF 

1F(PASS.EQ. 1)  GOTO  60 
C 

C - SELECT  FIN  TYPE - 

C 

IF(MODE.EQ. • 1* )  PRINT  1000 
lF(MODE.EQ. *2 • )  PRINT  1001 
50  PRINT  2200 
READ  9000, FIN 

IF(FIN.NE. '1* .AND. FIN. NE. •2* .AND. FIN. NE. '3’ )  THEN 
PRINT  1200 
GOTO  50 
ENDIF 

IF(PASS.EQ.1.AND.FIN.EQ. •!• )  GOTO  70 
IF(PASS.EQ.1.AND.FIN.EQ. •2* )  GOTO  75 
IF(PASS.EQ.1.AND.FIN.EQ. -a* )  GOTO  80 
C 

C -  SELECT  PROBLEM  TYPE  - 

C 

lF(MODE.EQ. • !• )  PRINT  1000 
IF<MODE.EQ. •2' )  PRINT  1001 
55  PRINT  2300 

READ  9000, TYPE 

IF(TYPE.NE. • 1' .AND.TYPE.NE. •2’ )  THEN 
PRINT  1200 
GOTO  55 
ENDIF 

IF(PASS.EQ. 1. AND. TYPE. EQ. • 1 • )  GOTO  98 
IF(PASS.EQ.1.AND.TYPE.EQ. •2' )  GOTO  96 


c- 

— 

-  INPUT  BASE 

TEMPERATURE 

c 

60 

IF(MODE.EQ. 

•1’) 

PRINT 

1000 

61 

IF(MODE.EQ. 
PRINT  3000 

•2'  ) 

PRINT 

1001 

READ  *,TB 

IF(UNITS.EQ. ■ . AND .TB .LT . -273 . 15D0 )  THEN 
PRINT  3001 
PRINT  1001 
GOTO  61 

ENDIF 

IF (UNITS. EQ. *2 • . AND . TB .LT . -460 . ODO )  THEN 
PRINT  3002 
PRINT  1001 
GOTO  61 
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ENDIF 

IF(PASS.EQ.l)  GOTO  105 
C 

C -  INPUT  FIN  LENGTH  - 

C 

65  IF(MODE.EQ. • 1* )  PRINT  1000 
lF(MODE.EQ. '2 • )  PRINT  1001 

66  PRINT  3100 
READ  *,L 

IF(L.LE.O.ODO)  THEN 
PRINT  3101 
PRINT  1001 
GOTO  66 
ENDIF 

IF(PASS.EQ.l)  GOTO  105 
C 

C -  INPUT  BASE  THICKNESS  IF  RECTANGLE  - 

C 

70  IF(FIN.EQ. • 1 • )  THEN 

IF(MODE.EQ. • 1* )  PRINT  1000 
IF(MODE.EQ. •2' )  PRINT  1001 

71  PRINT  3200 
READ  *,DEL 

IF(DEL.LE.O.ODO)  THEN 
PRINT  3201 
PRINT  1001 
GOTO  71 
ENDIF 
ENDIF 

IF  (PASS.EQ.l)  GOTO  105 
C 

C -  INPUT  BASE  AND  TIP  THICKNESSES  IF  TRAPEZOID 

C 

75  IF(FIN.EQ. '2' )  THEN 

IF(MODE.EQ. • 1 ’ )  PRINT  1000 
IF(MODE.EQ. '2' )  PRINT  1001 

76  PRINT  3300 

READ  *,DEL0,DELE 

IF ( DELO . LE . 0 . ODO , OR . DELE .LE . 0 . ODO )  THEN 
PRINT  3301 
PRINT  1001 
GOTO  76 
ENDIF 

IF ( DELO. LE. DELE)  THEN 
PRINT  3302 
PRINT  1001 
GOTO  76 
ENDIF 
ENDIF 

IF(PASS.EQ.l)  GOTO  105 
C 

c -  INPUT  BASE  THICKNESS  IF  TRIANGLE  - 

C 
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80  IF(FIN.EQ. 'S' )  THEN 

IF(MODE.EQ. • 1' )  PRINT  1000 
IF(MODE.EQ. -a* )  PRINT  1001 

81  PRINT  3200 
READ  *,DEL0 

IF(DELO.IiE.O.ODO)  THEN 
PRINT  3201 
PRINT  1001 
GOTO  81 
END  IF 
ENDIF 

IF(PASS.EQ. 1)  GOTO  105 
C 

C -  INPUT  CONDUCTIVITY  - 

C 

85  IF(MODE.EQ. • 1* )  PRINT  1000 
IF(MODE.EQ. *2 • )  PRINT  1001 

86  PRINT  3400 
READ  *,K 

IF(K.LE.O.ODO)  THEN 
PRINT  3401 
PRINT  1001 
GOTO  86 
ENDIF 

IF(PASS.EQ. 1)  GOTO  105 
C 

C -  INPUT  ABSORPTIVITY  - 

C 

90  1F(M0DE.EQ. • 1 • )  PRINT  1000 
IF(MODE.EQ. •2' )  PRINT  1001 

91  PRINT  3500 
READ  *,ABS 

IF(ABS.LT.O.ODO.OR.ABS.GT.l.ODO)  THEN 
PRINT  3501 
PRINT  1001 
GOTO  91 
ENDIF 

IF(PASS.EQ. 1)  GOTO  105 
C 

C -  INPUT  EMISSIVITY  - 

C 

92  IF(MODE.EQ. • 1 • )  PRINT  1000 
IF(MODE.EQ. '2' )  PRINT  1001 

93  PRINT  3600 
READ  *,EMIS 

IF(EMIS.LT.O.ODO.OR.EMIS.GT.l.ODO)  THEN 
PRINT  3601 
PRINT  1001 
GOTO  93 
ENDIF 

IF(PASS,EQ. 1)  GOTO  105 
C 

C - INPUT  EXTERNAL  HEAT  FLUX - 
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c 

94  IF(MODE.EQ. • 1* )  PRINT  1000 
IF(MODE.EQ. •2* )  PRINT  1001 

95  PRINT  3700 
READ  *,E 

IF(E.LT.O.ODO)  THEN 
PRINT  3701 
PRINT  1001 
GOTO  95 
END  IF 

IF(PASS.EQ. 1)  GOTO  105 
C 

C -  INPUT  HEAT  ENTERING  BASE  IF  SYNTHESIS  PROBLEM  - 

C 

96  IF(TYPE.EQ. •2* )  THEN 

IF(MODE.EQ. ' 1 • )  PRINT  1000 
IF(MODE.EQ. *2 • )  PRINT  1001 

97  PRINT  3800 
READ  *,Q 

IF(Q.LE.O.ODO)  THEN 
PRINT  3801 
PRINT  1001 
GOTO  97 
ENDIF 
END  IF 

IF(PASS.EQ.l)  GOTO  105 
C 

C -  INPUT  FIN  HEIGHT  IF  ANALYSIS  PROBLEM  - 

C 

98  IF(TYPE.EQ. • 1’ )  THEN 

IF(MODE.EO. • 1* )  PRINT  1000 
IF(MODE.EQ. *2 • )  PRINT  1001 

99  PRINT  3900 
READ  *,HT 

if(ht.le.o.odO)  then 

PRINT  3901 
PRINT  1001 
GOTO  99 
ENDIF 
ENDIF 

IF(PASS.EQ.l)  GOTO  105 
C 

CIE  *****  SUMMARIZE  DATA  INPUT  *************************************** 
C 

C -  PRINT  INPUTS  TO  SCREEN  - 

C 

105  IF(MODE.EQ. • 1 • )  PRINT  1000 
IF(MODE.EQ. •2* )  PRINT  1001 
PRINT  4000 
C 

C -  PRINT  UNITS  TYPE  - 

C 

IF(UNITS.EQ, • 1' )  PRINT  4005 
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IF(UNITS.EQ. *2' )  PRINT  4010 
C 

C - PRINT  FIN  TYPE - 

C 

IF(FIN.EQ. • !• )  PRINT  4015 
IF(FIN.EQ. •2' )  PRINT  4020 
IF(FIN.EQ. 'S* )  PRINT  4025 
C 

C -  PRINT  PROBLEM  TYPE  - 

C 

IF(TYPE.EQ. • 1* )  PRINT  4030 
IF(TYPE.EQ. '2 • )  PRINT  4035 
C 

C -  PRINT  COMMON  INPUTS:  BASE  TEMP  &  LENGTH  - 

C 

IF(UNITS.EO. • 1* )  PRINT  4040,  TB,L 
IF(UNITS.EQ. •2* )  PRINT  4041,  TB,L 
C 

C -  IF  RECTANGLE  OR  TRIANGLE  PRINT  BASE  THICKNESS  — 

C 

IF(FIN.EQ. • 1* .AND. UNITS. EQ. •!• )  PRINT  4045, DEL 
IF(FIN.EQ. • 1' .AND.UNITS.EQ. *2* )  PRINT  4046, DEL 
1F(FIN.EQ. •3' .AND.UNITS.EQ. ‘I* )  PRINT  4045, DELO 
IF(FIN.EQ. O’ .AND.UNITS.EQ. *2* )  PRINT  4046, DELO 
C 

C -  IF  TRAPEZOID  PRINT  BASE  &  TIP  THICKNESS  - 

C 

IF(FIN.EQ. •2' .AND.UNITS.EQ. •!• )  PRINT  4050 , DELO , DELE 
IF{FIN.EQ. ’2’ .AND.UNITS.EQ. ’2’)  PRINT  4051, DELO, DELE 


C 

C -  PRINT  COMMON  INPUTS:  CONDUCTIVITY,  ABSORPTIVITY 

C  EMISSIVITY,  AND  EXTERNAL  FLUX 

C 


IF(UNITS.EQ. ’ 1* )  PRINT  4055,K,ABS,EMIS,E 
IF{UNITS.EQ. ’2’ )  PRINT  4056,K,ABS,EMIS,E 
C 

C -  IF  SYNTHESIS  PROBLEM  PRINT  HEAT  INPUT  THRU  BASE 

C 

IF(TYPE.EQ. ’2 • .AND.UNITS.EQ. • 1’ )  PRINT  4060, Q 
IF(TYPE.EQ. ’2' .AND.UNITS.EQ. ’2' )  PRINT  4061, Q 
C 

-  IF  ANALYSIS  PROBLEM  PRINT  FIN  HEIGHT  - 

1F(TYPE.EQ. • 1 • .AND.UNITS.EQ. • 1’ )  PRINT  4065, HT 
IF(TYPE.EQ. • 1’ .AND.UNITS.EQ. ’2’ )  PRINT  4066, HT 
C 

C - CHECK  TO  SEE  IF  INPUTS  ARE  CORRECT - 

C 

110  PRINT  4700 

READ  9000, ANS 

1F(ANS.NE. • Y' .AND.ANS.NE. 'N’ )  THEN 
PRINT  1200 
GOTO  110 
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ENDIF 

IF(ANS.EQ. ‘y )  GOTO  115 


C 

C- 

C 


CORRECT  ANY  WRONG  INPUTS 


PASS  =  1 
PRINT  4750 
READ  9000, ANS 
IF(ANS.EQ. • !• ) 
IF(ANS.EQ. *2’ ) 
IF(ANS.EQ. 'S’ ) 
IF(ANS.EQ. 'A') 
IF(ANS.EQ. ’5* ) 


!•) 

GOTO 

44 

2’) 

GOTO 

50 

3-) 

GOTO 

55 

4-  ) 

GOTO 

60 

5') 

GOTO 

65 

6* .AND.FIN.EQ. •!• ) 

GOTO 

70 

6 • .AND.FIN.EQ. *2 • ) 

GOTO 

75 

6’  .AND.FIN.EQ.  O'  ) 

GOTO 

80 

T) 

GOTO 

85 

8') 

GOTO 

90 

9') 

GOTO 

92 

10'  ) 

GOTO 

94 

11* .AND.TYPE.EQ. '2* ) 

GOTO 

96 

11 • .AND.TYPE.EQ. 'I' ) 

GOTO 

98 

IF(ANS.EQ. ’7* ) 
IF(ANS.EQ. ’B' ) 
IF(ANS.EQ. 'S’ ) 
IF(ANS.EQ. • 10' ) 


PRINT  1200 
GOTO  105 
C 

ClF  *****  GOTO  SUBROUTINES  ******************************************* 
C 

c - for  ENGLISH;  CONVERT  ALL  TO  SI - 

C 

115  IF(UNITS.EO. '2' )  CALL  SI(TB,TE,L,HT, DEL, DELO, DELE, K,E,QI,Q) 

C 

-  CONVERT  C  TO  K  - 


C- 

C 

C 

C- 

C 


TB  =  TB  +  273.15D0 
— —  SELECT  SUBROUTINE 


IF(FIN.EQ. '1' .AND. TYPE. EQ. )  CALL  RTAN(TB,TE,L, HT,DEL,K, ABS, 
&  EMIS,E,QI,Q,EFF,FLAG) 

IF(FIN.EQ. ’ !• .AND.TYPE.EQ. •2' )  CALL  RTSY ( TB, TE, L, HT, DEL, K, ABS , 
&  EMIS,E,QI,Q,EFF,FLAG) 

IF(FIN.EQ. '2' .AND.TYPE.EQ.’l' )  CALL  TPAN(TB, TE,L, HT, DELO , DELE , 
&  K, ABS, EMIS,E,QI,Q,EFF, FLAG) 

IF(FIN.EQ. •2* .AND.TYPE.EQ. •2’ )  CALL  TPSY (TB, TE,L, HT, DELO , DELE , 
&  K, ABS, EMIS,E,QI,Q,EFF, FLAG) 

IF(FIN.EQ. '3' .AND.TYPE.EQ. •!• )  CALL  TRAN (TB, TE , L, HT, DELO , DELE , 
&  K, ABS, EMIS,E,OI,0,EFF, FLAG) 


IF(FIN.EQ. ’3' .AND.TYPE.EQ. •2’ )  CALL  TRSY ( TB, TE, L, HT, DELO , DELE, 
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&  K,ABS,EMIS,E,QI,Q,EFF,FLAG) 


C 

ClG  *****  SUMMARIZE  DATA  OUTPUT  ************************************** 


C 

C -  IF  THERE  IS  NO  SOLUTION,  START  OVER  WITH  NEW  PROBLEM 

C 

IF(FLAG.EQ. 1)  GOTO  125 
C 

C - CONVERT  K  TO  C - 

C 

TB  =  TB  -  273.15D0 
TE  =  TE  -  273.15D0 


C 


C - FOR  ENGLISH:  COVERT  ALL  BACK  TO  ENGLISH - 

C 

IF(UNITS.EQ. ’2’ )  CALL  ENG (TB,TE,L, HT, DEL, DELO , DELE , K, E, QI , Q) 
C 

c -  REPEAT  INPUTS  - 

C 


116  IF(MODE.EQ. • 1 • )  PRINT  1000 
IF(MODE.EQ. •2* )  PRINT  1001 
PRINT  4800 


C 

C -  PRINT  UNITS  TYPE  - 

C 

IF(UNITS.EQ. • 1* )  PRINT  4005 
IF(UNITS.EQ. •2' )  PRINT  4010 
C 

C - PRINT  FIN  TYPE - 

C 

IF(FIN.EQ. • !• )  PRINT  4015 
IF(FIN.EQ. *2 • )  PRINT  4020 
IF(FIN.EQ. *3' )  PRINT  4025 
C 

C -  PRINT  PROBLEM  TYPE  - 

C 


IF(TYPE.EQ. • 1' )  PRINT  4030 
IF(TYPE.EQ, ’2 ’ )  PRINT  4035 
C 

C -  PRINT  COMMON  INPUTS;  BASE  TEMP  &  LENGTH  - 

C 

IF(UNITS.EQ. • 1' )  PRINT  4040, TB,L 
IF(UNITS.EQ. •2' )  PRINT  4041, TB,L 
C 

C -  IF  RECTANGLE  OR  TRIANGLE  PRINT  BASE  THICKNESS 

C 


IF(FIN.EQ. • 1 • .AND. UNITS. EQ. • 1* )  PRINT  4045, DEL 
IF<FIN.EQ. • 1 • .AND.UNITS.EQ. '2 • )  PRINT  4046, DEL 
IF(FIN.EQ. •3' .AND. UNITS. EQ. • 1* )  PRINT  4045, DELO 
IF(FIN.EQ. •3* .AND.UNITS.EQ. ’2’ )  PRINT  4046, DELO 
C 

C - IF  TRAPEZOID  PRINT  BASE  &  TIP  THICKNESS - 

C 
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IF(FIN.EQ. *2 • .AND.UNITS.EQ. • !• )  PRINT  4050 , DELO , DELE 
IF(FIN.EQ. •2' .AND.UNITS.EQ. •2* )  PRINT  405 1 , DELO , DELE 


C 

C -  PRINT  COMMON  INPUTS;  CONDUCTIVITY,  ABSORPTIVITY 

C  EMISSIVITY,  AND  EXTERNAL  FLUX 

C 


IF(UNITS.EQ. • !• )  PRINT  4055 , K, ABS, EMIS , E 
IF(UNITS.EQ. •2‘ )  PRINT  4056 ,K, ABS, EMIS , E 
C 

C -  IF  SYNTHESIS  PROBLEM,  PRINT  HEAT  INPUT  - 

C 

IF(TYPE.5  sj. '2  •  .AND.UNITS.EO- •  1*  )  PRINT  4060, Q 
IF{TYPE.EQ. ’2' .AND.UNITS.EQ. •2* )  PRINT  4061, Q 
C 

C -  IF  ANALYSIS  PROBLEM,  PRINT  FIN  HEIGHT  - 

C 

IF(TYPE.EQ. ■ 1* .AND.UNITS.EQ. •!• )  PRINT  4065, HT 
IF(TYPE.EQ. .AND.UNITS.EQ. •2' )  PRINT  4066, HT 
C 


C -  SUMMARIZE  OUTPUT  - 

C 

PRINT  4805 
C 

C - IF  ANALYSIS  PROBr.EM,  PRINT  HEAT  OUTPUT 


C 

IF(TYPE.EQ. • 1' .AND.UNITS.EQ. •!• )  PRINT  4062, Q 
IF(TYPE.EQ. • 1 • .AND.UNITS.EQ. '2* )  PRINT  4063, Q 
IF(TYPE.EQ. ’ 1 • .AND.Q.LE.O.ODO)  PRINT  4064 
C 

C -  IF  SYNTHESIS  PROBLEM,  PRINT  FIN  HEIGHT  - 

C 

IF(TYPE.EQ.  ■2'  .AND.UNITS.EQ. ’I*  )  PRINT  406"', HT 
IF(TYPE.EQ. •2' .AND.UNITS.EQ. •2’ )  PRINT  4068, HT 
IF(TYPE.EQ. '2' .AND.HT.LE.O.ODO)  PRINT  4069 
C 

C -  PRINT  TIP  TEMPERATURE  - 

C 

IF(UNITS.EQ. • 1' )  PRINT  4810, TE 
IF(UNITS.EQ. '2* )  PRINT  4811,TE 

IF(UNITS.EQ. • 1 ■ .AND.TE.LT.-273. 15D0)  PRINT  4812 
IF(UNITS .EQ. • 2 • .AND.TE .LT.-460.0D0)  PRINT  4813 
C 

C -  PRINT  FIN  EFFICIENCY  - 

C 

PRINT  4815,EFF 

IF(EFF.LT.O.ODO.OR.EFF.GT.l.ODO)  PRINT  4816 
C 

c - ALSO  WRITE  OUTPUTS  TO  FILE - 

C 

WRITE (9, 4800) 

IF(UNITS.EQ. • 1' )  WRITE(9,4005) 

IF(UNITS.EQ. *2 ' )  WRITE(9,4010) 

IF(FIN.EQ. ■ 1 ’ )  WRITE(9,4015) 
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IF(FIN.EQ. '2 • ) 

IF{FIN.EQ. '3' ) 

IF(TyPE.EQ. • 1 • ) 

IF(TYPE.EQ. •2’ ) 

IF(UNITS.EQ. ■ 1' ) 

IF(UNITS.EQ. *2 • ) 

IF(FIN.EQ. 'I' .AND. UNITS. EQ. 
1F(FIN.EQ. '1' .AND. UNITS. EQ. 
1F(FIN.EQ. -3' .AND. UNITS. EQ. 
IF(FIN.EQ. 'a* .AND. UNITS. EQ. 
1F(FIN.EQ. '2 • .AND. UNITS. EQ. 
IF(FIN.EQ. *2 • .AND. UNITS. EQ. 
IF(UNITS.EQ. • 1* ) 

IF(UNITS.EQ. *2 • ) 

IF(TyPE.EQ. '2 • .AND. UNITS. EQ. 
IF(TyPE.EQ. '2 • .AND. UNITS. EQ. 
IF(TyPE.EQ. • 1 • .AND. UNITS. EQ. 
IF{TyPE.EQ. • 1* .AND. UNITS. EQ. 


) 

) 

) 

) 

) 

) 

WRITE (9 
WRITE ( 9 
-!•) 

'2*) 

■!•) 

•2* ) 


IF(TyPE.EQ. • 1* .AND. UNITS. EQ. *  1* ) 
IF(TyPE.EQ. • 1 • .AND. UNITS. EQ. •2’ ) 
IF(Q.LE.O.ODO) 

IF(TyPE.EQ. •2* .AND. UNITS. EQ. 'I' ) 
IF(TyPE.EQ. *2* .AND. UNITS. EQ. •2' ) 
IF(HT.LE.O.ODO) 

IF(UNITS.EQ. • 1 • ) 

IF(UN1TS.EQ. •2* ) 

IF ( UNITS. EQ. . AND . TE .LT. -273 . 15D0 ) 
IF ( UNITS. EQ. •2' .AND.TE.LT.-460.0D0) 


WRITE (9, 4020) 
WRITE(9,4025) 

WRITE (9, 4030) 

WRITE (9, 4035) 

WRITE (9, 4040)  TB,L 
WRITE{9,4041)  TB,L 
WRITE (9, 4045)  DEL 
WRITE (9, 4046)  DEL 
WRITE(9,4045)  DELO 
WRITE (9, 4 046)  DELO 
WRITE (9, 4050)  DELO, DELE 
WRITE(9,4051)  DELO, DELE 
,4055)  K,ABS,EMIS,E 
,4056)  K,ABS,EMIS,E 
WRITE(9,4060)  Q 
WRITE (9, 4061)  Q 
WRITE(9,4065)  HT 
WRITE (9, 4066)  HT 
WRITE (9, 4805) 

WRITE (9, 4062)  Q 
WRITE(9,4063)  Q 
WRITE (9, 4064) 

WRITE (9, 4067)  HT 
WRITE(9,4068)  HT 
WRITE (9, 4069) 
WRITE(9,4810)  TE 
WRITE(9,4811)  TE 
WRITE(9,4812) 
WRITE(9,4812) 

WRITE (9, 4815)  EFF 


-  CONVERT  OUTPUT  TO  OTHER  UNITS  ?  - 

PRINT  5199 
READ  9000, ANS 

IF{ANS.NE. 'y ' .AND.ANS.NE. ‘N’ )  THEN 
PRINT  1200 
GOTO  119 
END  IF 

IF(ANS.EQ. 'N' )  GOTO  120 
IF(UN1TS.EQ. ’ !• )  THEN 

CALL  ENG(TB,TE,L,HT, DEL, DELO, DELE, K,E,QI,Q) 
UNITS  =  '2' 

GOTO  116 
ENDIF 

IF(UN1TS.EQ. '2 • )  THEN 

CALL  SI (TB,TE,L,HT, DEL, DELO, DELE, K,E,QI,Q) 
UNITS  =  • 1' 

GOTO  116 
ENDIF 


DO  SAME  PROBLEM  AGAIN  WITH  DIFFERENT  INPUTS  7 


120 


c 

c - 

c 

125 


C 

ClH 

C 

800 

900 


1000 

1001 

1100 

1200 

1300 


1301 


1400 

1500 


PRINT  5200 
READ  9000, ANS 

IF(ANS.NE. 'Y' .AND.ANS.NE. ‘N* )  THEN 
PRINT  1200 
GOTO  120 
ENDIF 

IF(ANS.EQ. ’Y* )  GOTO  105 

-  NEW  PROBLEM  ?  - 

PRINT  5300 
READ  9000, ANS 

IF(ANS.NE. ’Y* .AND.ANS.NE. 'N* )  THEN 
PRINT  1200 
GOTO  125 
ENDIF 

IF(ANS.EQ. 'Y* )  GO  TO  44 
CLOSE (UNIT=9) 

STOP 

*****  FORMAT  STATEMENTS  ****************************************** 


FORMAT(/,'  ARE  YOU  IN  “CAPS  LOCK"  MODE  (Y  OR  N)  ?  •) 

FORMAT  (  /  ,  •  SELECT  VIEWING  MODE  { 1  OR  2  )  :  •  ,  /  , 

&  •  1.  SCREEN  VIEWING  (SCREEN  IS  CLEARED  BEFORE  NEXT  DISPLAY)  ',/, 
&  •  2.  LINE  VIEWING  (OUTPUT  TO  SCREEN  IS  LINE  BY  LINE)  ’) 

FORMAT(24(/) ) 

FORMAT ( / ) 

FORMAT (/,’  DO  YOU  WISH  TO  SKIP  THE  INTRODUCTION  AND  PROCEED'  , 

&  •  DIRECTLY  TO  A  PROBLEM',/, '  (Y  OR  N)  ?  ') 

FORMAT(/,'  ***************  INCORRECT  ENTRY!  ***************  ',/) 

FORMAT ( / , '  INTRODUCTION  ' , / , 

i 

&  '  THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  ANALYZE  AND  SYNTHESIZE  ',/, 

&  '  RADIATIVE  HEAT  TRANSFER  IN  LONGITUDINAL  FINS  OF  THREE  ' , / , 

&  '  PROFILES :  ' , / , 

«■  '  ',f, 

&  '  1 .  RECTANGULAR  ' , / , 

«.  •  2.  TRAPEZOIDAL  ',/, 

&  '  3 .  TRIANGULAR  ' , / , 

i 

&  •  TWO  TYPES  OF  PROBLEMS  ARE  SOLVED;  ',/, 

&  ■  '  , 

i  '  1.  ANALYSIS  ',/, 

5  '  2 .  SYNTHESIS  '  ) 

FORM^T(/,•  WTD  TWO  ENVIRONMENTS  ARE  CONSIDERED:  ',/, 

•  •  ,/, 

6  ■  1 .  FREE  SPACE  ' , / , 

Sr  '  2 .  NON-FREE  SPACE  '  ) 

FORMAT (/,'  -  PRESS  RETURN  TO  CONTINUE  -  ') 

FORMAT (/,'  FOR  THE  FIN  ANALYSIS  PROBLEM  THE  FOLLOWING  INPUTS’  ,/, 
&  '  ARE  NEEDE'1;  '  ,  /  , 

&  • 
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1600 


1700 


1800 


1900 


2000 


&  •  RECTANGLE  TRAPEZOID  TRIANGLE  • , / , 

i  .  -  -  - 

&  '  BASE  TEMPERATURE  BASE  TEMPERATURE  BASE  TEMPERATURE 

&  ■  LENGTH  LENGTH  LENGTH  • , / , 

&  '  BASE  THICKNESS  BASE  THICKNESS  BASE  THICKNESS  ' , / , 

&  •  TIP  THICKNESS  *  •,/, 

&  •  ABSORPTIVITY  ABSORPTIVITY  ABSORPTVITY  • , / , 

&  '  EMISSIVITY  EMISSIVITY  EMISSIVITY  *,/, 

&  •  CONDUCTIVITY  CONDUCTIVITY  CONDUCTIVITY 

&  *  HEIGHT  HEIGHT  HEIGHT  *,/, 

&  •  INCIDENT  FLUX  INCIDENT  FLUX  INCIDENT  FLUX  *,/, 

&  •  ' , 

&  •  *  TIP/BASE  THICKNESS  RATIO  ASSUMED  =0.01  ') 

FORMAT (/,*  THE  FIN  ANALYSIS  OUTPUTS  ARE:  *,/, 

&  '  'ill 

&  •  RECTANGLE  TRAPEZOID  TRIANGLE  • , / , 

4  .  -  -  - 

&  •  HEAT  DISSIPATED  HEAT  DISSIPATED  HEAT  DISSIPATED 

&  •  TIP  TEMPERATURE  TIP  TEMPERATURE  TIP  TEMPERATURE  • , / , 

&  •  EFFICIENCY  EFFICIENCY  EFFICIENCY  •  ) 

FORMAT ( / , •  FOR  THE  FIN  SYNTHESIS  PROBLEM  THE  FOLLOWING  INPUTS •  , / , 
&  •  ARE  NEEDED ;  • , / , 

«■  •  *,/, 

<1  •  RECTANGLE  TRAPEZOID  TRIANGLE  •  ,  / , 

4  .  -  -  - 

&  •  BASE  TEMPERATURE  BASE  TEMPERATURE  BASE  TEMPERATURE  ' , 

&  •  LENGTH  LENGTH  LENGTH  * , / , 

&  •  BASE  THICKNESS  BASE  THICKNESS  BASE  THICKNESS 

&  •  TIP  THICKNESS  * 

&  •  CONDUCTIVITY  CONDUCTIVITY  CONDUCTIVITY  • , / , 

&  •  ABSORPTIVITY  ABSORPTIVITY  ABSORPTIVITY 

&  •  EMISSIVITY  EMISSIVITY  EMISSIVITY 

&  •  HEAT  INPUT  HEAT  INPUT  HEAT  INPUT  * , / , 

i  •  'ill 

4  •  *  TIP/BASE  THICKNESS  ASSUMED  =0.01  •) 

FORMAT </,’  THE  FIN  SYNTHESIS  OUTPUTS  ARE:  *,/, 

•  '  ll  ! 

&  •  RECTANGLE  TRAPEZOID  TRIANGLE  *,/, 

4  .  -  -  -  .,/, 

&  •  HEIGHT  HEIGHT  HEIGHT 

4  ,  .pjp  temperature  TIP  TEMPERATURE  TIP  TEMPERATURE  •,/, 

&  •  EFFICIENCY  EFFICIENCY  EFFICIENCY  • ) 

FORMAT ( / , ■  THERE  ARE  TWO  MEASUREMENTS  MODES :  ’ , / , 

■  '  ih 

Jr  •  1.  SI  UNITS  •  ,/, 

i  '  2 .  ENGLISH  UNITS  •  ) 

FORMAT ( / , •  NOTES :  • , / , 

•  •  ,/, 

i  '  1.  MUST  BE  IN  "CAPS  LOCK"  MODE. 

i  •  2.  TO  ESCAPE,  PRESS  "CONTROL  BREAK"  AND  REENTER  PGM. 

Jr  '  3.  PROGRAM  INTERRUPTION  BY  "UNDERFLOW" , "OVERFLOW"  OR 

Jr  •  OTHER  PROCESSOR  ERRORS  INDICATES  PROBLEMS  IN  THE  ’  ,  / , 

Jr  •  INPUT  DATA  -  CHECK  FOR  CORRECTNESS .  '  ,  / , 
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2001 

2100 

2200 

2300 

3000 

3001 

3002 

3100 

3101 

3200 

3201 

3300 

3301 

3302 

3400 


&  •  4.  THE  MAXIMUM  NUMBER  OF  INTERVALS  SEARCHED  IS  100.  THIS 

&  •  CAN  BE  INCREASED  IN  INCREMENTS  OF  100  BY  MENU  CHOICE. 

&  •  THE  MAXIMUM  NUMBER  OF  ITERATIONS  FOR  FINDING  THE  ROOT  • , / , 

&  '  IS  50.  THIS  CAN  ALSO  BE  INCREASED  IN  INCREMENTS  OF 

it  •  50.  MOST  PROBLEMS  WILL  CONVERGE  WITHIN  THESE  LIMITS.  *,/, 

&  •  HOWEVER  SOME  INPUTS  CAUSE  DIVERGENCE  OR  OSCILLATION.  *,/, 


i  •  5.  ALTHOUGH  SIX  DECIMAL  PLACES  ARE  SHOWN  IN  THE  INPUT/  *,/, 

&  •  OUTPUT  DATA  SUMMARIES,  IT  IS  NOT  GUARANTEED.  IT 

&  •  DEPENDS  ON  THE  ACCURACY  OF  INPUT.  THE  BEST  THAT  *,/, 

&  '  CAN  BE  EXPECTED  IN  SINGLE  PRECISION  IS  5 .  '  ) 

FORMAT ( •  6 .  OUTPUT  IS  PRINTED  TO  THE  SCREEN  AND  TO  A  FILE  • , / , 

&  '  CALLED  "FIN.DAT".  UPON  EXITING  THE  PROGRAM  THIS 

&  •  CAN  BE  EDITED  WITH  YOUR  EDITOR  OR  WORD  PROCESSOR  * , / , 

&  *  OR  CAN  BE  PRINTED  DIRECTLY  WITH  THE  DOS  COMMAND 

&  •  "PRINT  FIN. DAT".  WARNING:  THIS  FILE  IS  WRITTEN 

&  •  OVER  WITH  EACH  BEGINNING  OF  THE  PROGRAM!  ’  ) 

FORMAT (/,'  SELECT  THE  APPROPRIATE  UNITS  (1  OR  2): 

4  • 

5  •  1.  SI  UNITS  (METERS, DEGREES  CENTIGRADE, WATT, ETC . )  •,/, 

6  •  2.  ENGLISH  UNITS  (FEET, DEGREES  FAHRENHEIT, BTU, ETC . )  •) 

FORMAT ( / , •  SELECT  THE  TYPE  OF  FIN  ( 1 ,  2 ,  OR  3 ) :  * , / , 

4  •  • ,/, 

it  •  1 .  RECTANGULAR  '  ,  / , 

&  •  2.  TRAPEZOIDAL  *,/, 

&  •  3 .  TRIANGULAR  *  ) 

FORMAT ( / , •  SELECT  THE  TYPE  OF  PROBLEM  ( 1  OR  2 ) :  * , / , 

4  •  •,/, 

i  •  1.  ANALYSIS  ’,/, 

4*2.  SYNTHESIS  ',/, 

4  •  •  ) 

FORMAT (•  WHAT  IS  THE  BASE  TEMPERATURE  (DEGREES  C  OR  F)  7  ’,/, 

4  •  EXAMPLE:  440  234.56  0.3504567E3  •) 

FORMAT (/,•  *»**  ERROR;  TEMPERATURE  MUST  BE  GREATER  THAN  OR’  , 

&  •  EQUAL  TO  ABSOLUTE’  ,/, 

&  ’  ZERO  (-273.15  DEGREES  C)  I  ’) 

FORMAT (/,’  ****  ERROR;  TEMPERATURE  MUST  BE  GREATER  THAN  OR’  , 

&  •  EQUAL  TO  ABSOLUTE ’  , / , 

&  •  ZERO  (-460  DEGREES  F)  !  ’) 

FORMAT (’  WHAT  IS  THE  LENGTH  OF  THE  FIN  (M  OR  FT)  ?  ’,/, 

&  ’  EXAMPLE:  1  3.4565  0.73456E1  ’) 

FORMAT (/,’  ****  ERROR;  LENGTH  MUST  BE  GREATER  THAN  ZERO  !  ’) 

FORMAT (’  WHAT  IS  THE  FIN  THICKNESS  AT  THE  BASE  (M  OR  FT)  7  ’,/, 

&  ’  EXAMPLE;  .0095  0.0034  0.7345645-2  ’) 

FORMAT (/,'  **•*  ERROR;  BASE  THICKNESS  MUST  BE  GREATER  THAN’  , 

&  ’  ZERO  !  ’  ) 

FORMAT ( ’  WHAT  IS  THE  FIN  THICKNESS  AT  THE  BASE  AND  TIP  7 ’  , 

&  ’  (M  OR  FT)  7’,/,’  ENTER  2  NUMBERS  SEPARATED  BY  A  SPACE.’  ,/, 

4  ’  EXAMPLE:  .034  0.0068  OR  0.34E-2  0.234567E-2  ’) 

FORMAT(/,’  ****  ERROR;  THICKNESS  MUST  BE  GREATER  THAN  ZERO  !  ’) 

FORMAT (/,'  ****  ERROR:  BASE  THICKNESS  MUST  BE  GREATER  THAN’  , 

&  ’  TIP  THICKNESS  ! ’  ) 

FORMAT(’  WHAT  IS  THE  CONDUCTIVITY  (W/M-C  OR  BTU/FT-HR-F)  7  ’,/, 

4  ’  EXAMPLE:  154  33.06  0.2563455E+3  •) 
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3401 

FORMAT ( / 

,  •  ****  ERROR:  CONDUCTIVITY  MUST  BE  GREATER  ZERO  1 

•) 

3500 

FORMAT ( • 

WHAT 

IS  THE  ABSORPTIVITY 

(DIMENSIONLESS)  7* 

&  '  EXAMPLE:  . 

12  0.56  2.5E-1 

•) 

3501 

FORMAT ( / 

,  •  ****  ERROR:  ABSORPTIVITY  MUST  BE  GREATER  THAN  OR* 

r 

&  ‘  EQUAL 

TO  ZERO* 

& 

1 

AND  LESS  THAN  OR  EQUAL  TO  ONE  I 

•) 

3600 

FORMAT ( • 

WHAT 

IS  THE  EMISSIVITY  (DIMENSIONLESS)  7 

&  '  EXAMPLE:  . 

4  0.55  75.0E-2 

•) 

3601 

FORMAT (/ 

,  •  ****  ERROR:  EMISSIVITY 

MUST  BE  GREATER  THAN  OR’ 

f 

&  '  EQUAL 

TO  ZERO* 

& 

• 

AND  LESS  THAN  OR  EQUAL  TO  ONE  1 

') 

3700 

FORMAT ( • 

WHAT 

IS  THE  EXTERNAL  HEAT 

FLUX  (W  OR  BTU/HR)  7 

&  •  EXAMPLE:  675  450.78  0.234567+3 

•) 

3701 

FORMAT ( / 

,  •  ****  ERROR:  EXTERNAL  BEAT  FLUX  MUST  BE  GREATER’ 

/ 

&  '  THAN 

OR  EQUAL  TO  ZERO  I  * 

) 

3800 

FORMAT ( • 

WHAT 

IS  THE  HEAT  ENTERING 

THE  BASE  (W  OR  BTU/HR)  7 

&  •  EXAMPLE:  555  204.56  0.750567+3 

•) 

3801 

FORMAT (/ 

,  •  ****  ERROR:  HEAT  INPUT 

THROUGH  THE  BASE  MUST  BE* 

F 

&  •  GREATER  THAN  ZERO  ! ’ 

) 

3900 

FORMAT ( • 

WHAT 

IS  THE  HEIGHT  OF  THE 

FIN  (M  OR  FT)  7 

’ f 

&  ‘  EXAMPLE:  . 

45  0.95  7.5E-1 

' ) 

3901 

FORMAT ( / 

,  •  ****  ERROR:  FIN  HEIGHT 

MUST  BE  GREATER  THAN  ZERO 

1  ') 

4000 

FORMAT ( / 
£.  • 

,  *  INPUT  DATA  SUMMARY: 

•) 

•) 

4005 

FORMAT ( • 

(1) 

UNITS 

=  SI 

4010 

FORMAT ( • 

(1) 

UNITS 

=  ENGLISH 

• ) 

4015 

FORMAT ( • 

(2) 

FIN 

=  RECTANGLE 

•) 

4020 

FORMAT ( • 

(2) 

FIN 

=  TRAPEZOID 

• ) 

4025 

FORMAT ( • 

(2) 

FIN 

=  TRIANGLE 

•) 

4030 

FORMAT ( • 

<3) 

PROBLEM 

=  ANALYSIS 

') 

4035 

FORMAT ( * 

(3) 

PROBLEM 

=  SYNTHESIS 

•) 

4040 

FORMAT ( • 

<4) 

BASE  TEMPERATURE 

=*,G15.6,*  C 

& 

(5) 

FIN  LENGTH 

=  * ,G15.6,  *  M 

•) 

4045 

FORMAT ( • 

<6) 

BASE  THICKNESS 

=* ,G15.6, *  M 

•  ) 

4050 

FORMAT ( • 

(6) 

BASE/TIP  THICKNESSES 

=  •  ,G15.6,2X,G15.6, *  M 

•  ) 

4055 

FORMAT ( • 

(7) 

CONDUCTIVITY 

=  * ,G15.6,  *  W/M-C 

& 

(8) 

ABSORPTIVITY 

=* ,G15.6, 

/, 

& 

(9) 

EMISSIVITY 

=* ,G15.6, 

!, 

& 

(10) 

EXTERNAL  HEAT  FLUX 

=* ,G15.6, *  W 

•  ) 

4060 

FORMAT ( • 

(11) 

HEAT  INPUT  THRU  BASE 

=* ,G15.6, *  W 

• ) 

4062 

FORMAT ( • 

(12) 

HEAT  INPUT  THRU  BASE 

=* ,G15.6, *  W 

•  ) 

4065 

FORMAT ( • 

(11) 

FIN  HEIGHT 

=* ,G15.6, *  M 

•) 

4067 

FORMAT ( • 

(12) 

FIN  HEIGHT 

=*,G15.6,*  M 

•  ) 

4041 

FORMAT ( • 

(4) 

BASE  TEMPERATURE 

=*,G15.6, *  F 

& 

(5) 

FIN  LENGTH 

=*,G15.6,*  FT 

' ) 

4046 

FORMAT ( • 

(6) 

BASE  THICKNESS 

=*,G15.6,*  FT 

• ) 

4051 

FORMAT ( ’ 

(6) 

BASE/TIP  THICKNESSES 

=  *  ,G15.6,2X,G15.6, *  FT 

• ) 

4056 

FORMAT ( • 

(7) 

CONDUCTIVITY 

=*,G15.6,*  BTU/HR-FT-F 

& 

(8) 

ABSORPTIVITY 

= ’ ,G15 . 6 , 

/, 

& 

(9) 

EMISSIVITY 

=• ,G15.6, 

/, 

& 

(10) 

EXTERNAL  HEAT  INPUT 

**,G15.6,*  BTU/HR 

'  ) 

4061 

FORMAT ( • 

(11) 

HEAT  INPUT  THRU  BASE 

**,G15.6,*  BTU/HR 

• ) 

4063 

FORMAT ( • 

(12) 

HEAT  INPUT  THRU  BASE 

=*,G15.6,*  BTU/HR 

•) 
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4064  FORMAT (•  ****  WARNING:  HEAT  INPUT  THROUGH  BASE  LESS  <0,*  , 

&  •  RECHECK  YOUR  INPUTS  I '  ) 

4066  FORMATC  (11)  FIN  HEIGHT  =*,G15.6,*  FT  •) 

4068  FORMATC  (12)  FIN  HEIGHT  =',G15.6,’  FT  •) 

4069  FORMAT( •  **»*  WARNING:  FIN  HEIGHT  <  OR  =  0,  RECHECK  YOUR*  , 

&  •  INPUTS  1 •  ) 

4700  FORMAT(/,'  IS  THIS  CORRECT  (Y  OR  N)  7  ') 

4750  FORMAT (/,'  WHAT  IS  THE  NUMBER  OF  THE  INCORRECT  ENTRY  7  ’) 

4800  FORMAT(/*  INPUT/OUTPUT  SUMMARY: 

6  • 

5  •  **********************  INPUT  ***************************  ■/) 

4805  FORMAT (/  , 

f,  >  **********************  OUTPUT  **************************  •/) 

4810  FORMATC  (13)  TIP  TEMPERATURE  =‘,G15.6,’  C  •) 

4811  FORMATC  (13)  TIP  TEMPERATURE  =*,G15.6,'  F  •) 

4812  FORMATC  ****  WARNING;  FIN  TIP  TEMP  <  -273.15  DEG  C,  RECHECK*  , 

6  *  YOUR  INPUTS  I  *  ) 

4813  FORMATC  ****  WARNING;  FIN  TIP  TEMP  <  -460.0  DEG  F,  RECHECK*  , 

&  *  YOUR  INPUTS  1  *  ) 

4815  FORMAT(*  (14)  FIN  EFFICIENCY  =*,G15.6  ) 

4816  FORMAT(*  ****  WARNING:  FIN  EFFICIENCY  <  0  OR  >  1,  RECHECK*  , 

&  *  YOUR  INPUTS  I  *  ) 

5199  FORMAT(/*  WOULD  YOU  LIKE  TO  SEE  OUPUT  IN  OTHER  UNITS  (Y  OR  N)  7*  ) 

5200  FORMAT(/*  WOULD  YOU  LIKE  TO  DO  SAME  PROBLEM  AGAIN  WITH*  , 

&  *  DIFFERENT  INPUTS  (Y  OR  N)  7*  ) 

5300  FORMAT(/*  WOULD  YOU  LIKE  TO  START  OVER  WITH  A  NEW  PROBLEM*  , 

&  *  (Y  OR  N)7*  ) 

9000  FORMAT(A2) 

END 

C 

C 


C2A  *****  RTAN:  HEADERl  ********************************************** 
C 

C  SUBROUTINE  RTAN ( TB, TE, HT, L, DEL, K,ABS, EMIS , E, QI , Q, EFF, FLAG) 

C 

C  FIN  PROFILE;  RECTANGLE 
C 

C  PROBLEM  TYPE;  ANALYSIS 
C 

C  INPUT:  TB,HT,L,DEL,K,ABS,EMIS,E 
C 

C  OUTPUT:  Q,TE,QI,EFF 
C 

c - 

c 

C  PARAMETERS : 

C 

C  TB  -  TEMPERATURE  AT  THE  BASE  OF  THE  FIN 
C  TE  -  TEMPERATURE  AT  THE  TIP  OF  THE  FIN 
C  HT  -  HEIGHT  OF  THE  FIN 
CL-  LENGTH  OF  THE  FIN 
C  DEL  -  THICKNESS  FIN 
C  K  -  CONDUCTIVITY  OF  THE  FIN 
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C  ABS  -  ABSORPTIVITY  OF  THE  FIN 

C  EMIS  -  EMISSIVITY  OF  THE  FIN 

C  E  -  EXTERNAL  HEAT  INCIDENT  ON  THE  FIN 

C  Q  -  HEAT  DISSIPATED  BY  THE  FIN 

C  EFF  -  EFFICIENCY  OF  THE  FIN 

C  FLAG  -  0  =  CONVERGENCE,  1  =  DIVERGENCE 

C 

C - 

C 

C  FUNCTIONS : 

C 

C  FCN1(X,TB,TE,HT,DEL,K,K1,K2) : 

C 

C  FUNCTION  WHOSE  ROOT  IS  FOUND  (BASE  DERIVATIVE  &  TIP  TEMP) 

C 

C - 

C 

C  SUBROUTINES : 

C 

C  MDLIN1(FCN1,X1,X2,XR,TB,HT,DEL,K,K1,K2) : 

C 

C  FINDS  THE  ROOT  OF  FUNCTION  FCNl  BY  THE  METHOD  OF  MODIFIED 
C  LINEAR  INTERPOLATION 
C 

C  RKFSY1(DERIV1,XO,TO,TEND,F,H,DEL,K,K1,K2,TOL) : 

c 

C  SOLVES  A  SECOND  ORDER  NONLINEAR  DE  BY  RUNGE-KUTTE-FELHBERG 
C  METHOD 
C 

DERIV1(T,F,DEL,K,K1,K2) 

C  COMPUTES  THE  DERIVATIVES  FOR  RKFSYl 
C 

C2B  *****  RTAN:  MAINl  ************************************************ 
C 

SUBROUTINE  RTAN (TB, TE, L, HT, DEL, K, ABS, EMIS, E, QI , Q, EFF, FLAG) 

REAL*8  TB,TE,L,HT,DEL,K,ABS,EMIS,E,QI,Q,EFF,X1,X2,XR, 

(  SB,Kl,K2,FCNl,Fl,F2,TBDER 

INTEGER  I, J, PASS, FLAG 
CHARACTER* 2  ANS 
EXTERNAL  FCNl 


C 

C -  DEFINE  CONSTANTS 

C 

SB  =  5.67051D-12 
K1  =  2,0D0*SB*EMIS 
R2  =  E*ABS 
C 

C -  SCALE  INPUTS  - 

C 

HT  =  HT*100.0D0 
DEL  =  DEL* 100. ODD 
L  *  L*100.0D0 
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K  =  K*0.01D0 
C 

C -  START  COUNTERS  - 

C 

PASS  =  0 
1  =  0 
J  =  1 
C 

c -  START  SEARCH  WITH  INTERVAL  [-1,0]  - 

C 

PRINT  200 
Xl  =  O.ODO 
X2  =  -l.ODO 

Fl  =  FCNl(Xl,TB,TE,HT,DEL,K,Kl,K2) 

5  1  =  1  +  1 
C 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  100  INTERVALS 

C 

IF(I.GT. J»100)  GOTO  10 
IF(PASS.EQ.O)  PRINT  100,1 
IF(PASS.EQ. 1)  PRINT  200 
IF(PASS.GE.l)  PRINT  300, X1,F1,X2,F2 


F2  =  FCN1(X2,TB,TE,HT,DEL,K,K1,K2) 

C 

C -  CHECK  FOR  FUNCTION  VALUES  OPPOSITE  IN  SIGN 

C 


IF (Fl*F2.GE. O.ODO)  THEN 
XI  =  X2 
Fl  =  F2 

X2  =  X2  -  l.ODO 
1F(PASS.GE. 1)  PASS  =  PASS  +  1 
GOTO  5 
ENDIF 
C 

C - IF  INTERVAL  FOUND,  GO  FIND  ROOT - 

C 

PRINT  201 
GOTO  20 
C 

C - IF  NO  INTERVAL  FOUND,  CONTINUE  SEARCH  ? - 

C 

10  PRINT  400,J*100 
READ  500, ANS 

IF(ANS.NE. 'Y* .AND.ANS.NE. ’N* )  THEN 
PRINT  600 
GOTO  10 
ENDIF 

IF(ANS.EQ. ’N' )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 
C 

C -  IF  UNABLE  TO  FIND  INTERVAL,  LOOK  AT  FUNCTION  VALUES  7 
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c 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. 'Y' .ANO.ANS.NE. ’M< )  THEN 
PRINT  600 
GOTO  15 
ENDIF 

IF(ANS.EQ. ‘N* )  THEN 
FLAG  =  1 
RETURN 
ENDIF 
PASS  s  1 
1  =  0 
J  =  1 

XI  =  O.ODO 
X2  =  -l.ODO 

FI  =  FCNl(Xl,TB,TE,HT,DEL,K,Kl,K2) 

GOTO  5 
C 

C -  GET  ROOT  IN  INTERVAL:  DERIV  AT  BASE  AND  TEMP  AT  TIP  - 

C 

20  CALL  MDLIN1(FCN1,X1,X2,XR,TB,TE,DEL,K,K1,K2,HT) 

TBDER  =  XR 
C 

C -  CALCULATE  IDEAL  HEAT  DISSIPATED  - 

C 

QI  =  2.0D0*SB*EMIS*HT*L*TB**4.ODO 
C 

C -  CALCULATE  REAL  HEAT  DISSIPATED  - 

C 

Q  =  -K*DEL*L*TBDER 
C 

C -  CALCULATE  FIN  EFFICIENCY  - 

C 

EFF  =  Q/QI 
C 

C -  SCALE  OUTPUTS  - 

C 

HT  =  HT*0.01D0 
DEL  =  DEL* 0.0 IDO 
L  =  L*0.01D0 
K  =  K*100.0D0 
RETURN 
C 

C -  FORMAT  - 

C 

100  FORMAT (•+',• LOOKING  IN  INTERVAL:  ’,13, •  •) 

200  FORMAT (/,'  COMMENCING  SEARCH  FOR  INTERVAL  TO  BRACKET  ROOT _  •,/) 

201  FORMAT (/,'  INTERVAL  FOUND,  WILL  NOW  LOOK  FOR  ROOT  _  ’,/) 

300  FORMATC  XL  = ' , 2X, DIO . 4 , 2X, 'F(XL)  »• ,2X,D10.4,2X, 

&  ’XR  =• ,2X,D10.4,2X, •F{XR)  =',2X,D10.4) 

400  FORMAT (/,'  NO  ROOT  FOUND  IN  ’,13, *  INTERVALS.  CONTINUE  SEARCH'  , 
4  •  (Y  OR  N)  ?•  ) 
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500  FORMAT (A2) 

600  FORMAT(/,'  ****************  INCORRECT  ENTRY  I  *****************  ') 

700  FORMAT</,'  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N)7 

&  •  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  IF 
&  '  F(XR)  IS  OF  THE  SAME  SIGN  AS  F(XL)  AND  ITS  ABSOLUTE  VALUE 
&  •  IS  INCREASING,  THERE  IS  LITTLE  PROBABILITY  OF  FINDING  AN  *,/, 
&  •  INTERVAL.  THE  SAME  GOES  FOR  VALUES  OF  F(XR)  THAT  OSCILLATE*,/, 
&  *  AROUND  A  CERTAIN  VALUE.  FOR  F(XR)  TO  BECOME  OPPOSITE  IN  *,/, 
&  *  SIGN  TO  F{XR),  IT  MUST  PASS  THROUGH  ZERO.  •) 

END 

C 

C2C******  RTAN:  FCNl (X,TB,TE,HT,DEL,K,Kl,K2)  ************************ 

C 

C  PARAMETERS : 

C 

C  X  -  INDEPENDENT  VARIABLE  (BASE  DERIV) 


c 

TB 

-  BASE  TEMP 

c 

TE 

-  TIP  TEMP 

c 

HT 

-  FIN  HEIGHT 

c 

DEL 

-  FIN  THICKNESS 

c 

K  - 

FIN  CONDUCTIVITY 

c 

Kl 

-  CONSTANT  =  2.0D0*SB*EMIS 

c 

K2 

-  CONSTANT  =  DABS*E 

C 

C - 

C 

REAL*8  FUNCTION  FCNl ( X, TB, TE,HT, DEL, K, K1 , K2 ) 

REAL*8  X,X0,XFINAL,T0(2) ,TEND(2) ,F(2) ,H,TB,TE,HT,DEL, 

&  K,Kl,K2,TSAVE,TOL 
EXTERNAL  DERIV 1 
C 

C -  INITIALIZE  LEFT  END  BC*S:  FCT  VALUE  KNOWN,  GUESS  DERIV  - 

C 

T0(1)=  TB 
T0(2)=  X 
C 

c - INITIALIZE  RIGHT  END  BC*S  TO  0 - 

C 

TEND( 1)=0.0D0 
TEND(2)=0.0D0 
C 

C -  DEFINE  PARAMETERS  - 

C 

TSAVE  =  1000. ODD 
H  =  O.lDO 
TOL  =  O.OOOlDO 
XO  =  O.ODO 
XFINAL  =  HT 

10  IF  (XO  .LT.  XFINAL+O.OOOIDO)  THEN 

CALL  RKFS Y 1 ( DERIV 1 , X  0 , T  0 , TEND , F , H , DEL , K , K 1 , K2 , TOL ) 

INCREASING  T  MEANS  DIVERGENCE  - 


C 

C 

c 
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IF  (TEND(l) .GE.TSAVE)  GO  TO  20 
T0(1)  =  TEND{1) 

TSAVE  =  TEND(l) 

TO (2)  =  TEND(2) 

GOTO  10 
ENDIF 

20  CONTINUE 

FCNl  =  TEND{2) 

TE  =  TEND(l) 

RETURN 

END 

C 

C2D  *****  RTAN;  RKFSyl(DERIVl,XO,TO,TEND,F,H,DEL,K,Kl,K2,TOL) 

C 

C  PARAMTERS : 

C 

C  RKFSYl  -  SUBROUTINE  THAT  SOLVES  A  SYSTEM  OF  2  FIRST  ORDER 
C  DIFFERENTIAL  EQUATIONS  BY  THE  RUNGE-KUTTA-FEHLBERG 

C  METHOD.  THE  EQUATIONS  ARE  OF  THE  FORM: 

C 

C  DT/DX  =  Y  =  F1(X,T) 

C  DY/DX  =  F2(X,T,Y) 

C 

C  DERIVl  -  A  SUBROUTINE  THAT  COMPUTES  VALUES  OF  THE  2  DERIVATIVES. 
C 

C  XO  -  THE  INITIAL  VALUE  OF  THE  INDEPENDENT  VARIABLE 
C  TO  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 
C  TEND  -  AN  ARRAY  THAT  RETURNS  THE  PINAL  VALUES  OF  THE  FUNCTIONS 
C  F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 

C  H  -  THE  INCREMENT  TO  T,  THE  STEP  SIZE 

C  DEL  -  FIN  THICKNESS 

C  K  -  FIN  CONDUCTIVITY 

C  K1  -  CONSTANT  =  2.0*SB*EMIS 
C  K2  -  CONSTANT  =  E*ABS 
C  TOL  -  TOLERANCE 

C  TWRK  -  AN  ARRAY  USED  TO  HOLD  INTERMEDIATE  VALUES  DURING  THE 
C  COMPUTATION.  IT  MUST  BE  DIMENSIONED  OF  SIZE  6X2. 

C 

C - 

c 

SUBROUTINE  RKFS Y 1 ( DERIVl , X 0 , TO , TEND , F , H , DEL , K , K 1 , K2 , TOL ) 
REAL*8  X0,T0(2) ,TEND(2) ,F(2),TWRK(6,2) ,H,K,K1,K2,T0L,DEL, 


&  ERROR, SUM, STOREH,XEND 

INTEGER  I 
C 

C - INITIALIZE  FOR  INTERVAL  [XO,  XO+H] 

C 


XEND  =  XO+H 
STOREH  *  H 
C 

C - CHECK  TO  SEE  IF  FINISHED 

C 

1  IF(XO.GE.XEND)  THEN 
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noo  ooo  non 


H  =  STOREH 
RETURN 
ELSE 
C 

C - GET  FIRST  ESTIMATE  OF  THE  DELTA  X*S - 

C 

CALL  DERIV1(T0,F,DEL,K,K1,K2) 

DO  10  I  =  1,2 

TWRK(1,I)  =  H*F(I) 

TEND(I)  =  T0(I)+TWRK(1,I)/4.0D0 
10  CONTI  VIE 

-  GET  SECOND  ESTIMATE  - 

CALL  DERIVl (TEND, F, DEL, K,K1,K2) 

DO  20  I  =  1,2 

TWRK(2,I)  =  H*F(I) 

TEND(I)  =  TO(I)+(TWRK(1,1)*3.0DO+TWRK(2,I)*9.0DO)/32.0DO 
20  CONTINUE 

-  GET  THIRD  ESTIMATE  - 

CALL  DERIV1(TEND,F,DEL,K,K1,K2) 

DO  30  I  =  1,2 

TWRK(3,I)  =  H*F(I) 

TEND<I)  =  T0(I)+(TWRK(1,1)*1932.0D0-TWRK(2,I)*7200.0D0 
&  +  TWRK(3,I)*7296.0D0)/2197.0D0 

30  CONTINUE 

-  get  fourth  ESTIMATE  - 

CALL  DERIV1(TEND,F,DEL,K,K1,K2) 

DO  40  I  =  1,2 

TWRK(4,I)  =  H*F(I) 

TEND(I)  =  T0(1)+(439.0D0*TWRK(1,I)/216.0D0-8.0D0»TWRK(2,I) 
&  +  3680.0D0*TWRK(3,I)/513.0D0-845.0D0*TWRK(4,1)/4104.0D0) 

40  CONTINUE 
C 

C - GET  FIFTH  ESTIMATE - 

C 

CALL  DERIVl (TEND, F, DEL, K,K1,K2) 

DO  50  I  =  1,2 

TWRK(5,I)  =  H*F(I) 

TEND(I)  =  T0(I)-8.0D0*TWRK(l,I)/27.0+2.0D0*TWRK(2,I) 


&  -  3544.0DO*TWRK(3,I)/2565.0DO+1859.0DO*TWRK(4,I)/4104.0DO 

&  -11.0D0*TWRK(5,I)/40.0D0 

50  CONTINUE 
C 

C -  GET  SIXTH  ESTIMATE  - 

C 


CALL  DERIVl ( TEND, F, DEL, K,K1,K2) 
DO  60  I  =  1,2 

TWRK(6,I)  =  H*F(I) 
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60  CONTINUE 


C 

C -  ESTIMATE  THE  ERROR  BY  COMPUTING  THE  DIFFERENCE  BETWEEN 

C  THE  FOURTH  AND  FIFTH  ORDER  EQUATIONS 

C 

SUM  =  O.ODO 
ERROR  =  O.ODO 
DO  70  I  =  1,2 

SUM  =  DABS(TWRK(1,I)/360.0D0-128.0D0*TWR.I(3,I)/4275.0D0 

&  -  2197.0DO*TWRK(4,I)/75240.0DO+TWRK(5,I)/50.0DO 

&  +  2 .0D0*TWRK(6,I) /55.0DO) 

ERROR  =:  DMAXl  (ERROR,  SUM) 

70  CONTINUE 
C 

C -  IF  ERROR  IS  LESS  THAN  TOLERANCE,  COMPUTE  X  AT  THE  - 

C  END  OF  THE  INTERVAL  FROM  A  WEIGHTED  AVERAGE  OF  THE  SIX 

C  ESTIMATES  AND  RETURN 

C 

IF ( ERROR. LT.TOL)  THEN 


DO  80  I  =  1,2 


TEND(I)  =  T0(I)+16.0D0*TWRK{1,I)/135.0D0 
&  +  6656.0D0*TWRK(3,I)/12825.0D0 

&  +  28561. 0D0*TWRK{4,I)/5643O. 0-9. OD0*TWRK(5,I)/5O.ODO 

&  +  2.0D0*TWRK(6,I)/55.0D0 

T0<I)  =  TEND(I) 

80  CONTINUE 
XO  =  XO+H 
ENDIF 
C 

c -  IF  ERROR  IS  GREATER  THAN  TOLERANCE,  DECREASE  STEP  &  - 

C  GO  AGAIN 

C 

IF < ERROR. GT.TOL)  THEN 
H  =  H/2.0D0 
ENDIF 
C 

C -  IF  ERROR  IS  SIGNIFICANTLY  LESS  THAN  TOLERANCE,  RELAX  - 

C  STEP 

C 

IF(ERROR  .LT.  H*TOL/ 10 . ODO )  THEN 
H  =  H*2.0D0 
ENDIF 
C 

C - IF  OVERSHOOT  END,  REDUCE  STEP - 

C 

IF(X0  +  H  .GT.  XEND)  THEN 
H  =  XEND-XO 
ENDIF 
C 

ENDIF 

C 

C 


GO  TO  1 
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END 


C 

C2E******  RTAN:  DERIVl(T,F,DEL,K,Kl,K2)  ****************************** 
C 

C  PARAMETERS : 

C 

C  T  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FtTOCTIONS 

C  F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 

C  DEL  -  FIN  THICKNESS 

C  K  -  FIN  CONDUCTIVITY 

C  K1  -  CONSTANT  =  2.0*SB*EMIS 

C  K2  -  CONSTANT  =  E*ABS 

C 

C - 

c 

SUBROUTINE  DERIVl ( T , F , DEL , K , Kl , K2 ) 

REAL*8  T(2) ,F(2) ,K,K1,K2,DEL 
F(l)  =  T(2) 

F(2)  =  K1*T(1)**4.0D0/(K*DEL)  -  K2/(K*DEL) 

RETURN 

END 

C 

C2F  *****  RTAN;  MDLINl (FCNl , XI , X2 , XR,TB,TE, DEL, K, Kl ,K2 , HT)  ********** 

• , /, ,TB,TE, DEL, ******* 

C 

C  PARAMETERS ; 

C 

C  FCNl  -  FUNCTION  THAT  COMPUTES  VALUES  FOR  F.  MUST  BE  DECLARED 
C  EXTERNAL  IN  CALLING  PROGRAM 

C  Xl,X2  -  INITIAL  VALUES  OF  X.  F(X)  MUST  CHANGE  SIGNS  AT  THESE 
C  POINTS 

C  XR  -  RETURNS  THE  ROOT  TO  THE  MAIN  PROGRAM 
C  TB  -  BASE  TEMPERATURE 
C  TE  -  TIP  TEMPERATURE 
C  DEL  -  FIN  THICKNESS 
C  K  -  FIN  CONDUCTIVITY 
C  Kl  -  CONSTANT  =  2.0*SB*EMIS 
C  K2  -  CONSTANT  =  E*ABS 
C  HT  -  FIN  HEIGHT 
C 

c - 

c 

SUBROUTINE  MDLINl ( FCNl , Xl , X2 , XR,TB, TE, DEL, K, Kl , K2 , HT ) 

REAL* 8  XERR,FSAVE,Fl,F2,FR,Xl,X2,XR,XTOL,FTOL,TB,TE,DEL, 

&  K,K1,K2,FCN1,HT 

INTEGER  I, J, PASS 
CHARACTER* 2  ANS 

DATA  XTOL,FTOL/0 . OOOlDO, 0 . OOOOlDO/ 

C 

C -  INITIALIZE  VALUES  - 

C 

PRINT  200 
XibAV  =  XI 
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X2SAV  =  X2 

FI  =  FCN1(X1,TB,TE,HT,DEL,K,K1,K2) 

F2  =  FCN1(X2,TB,TE,HT,DEL,K,K1,K2) 

FlSAV  =  Fl 
F2SAV  =  F2 
FSAVE  =  Fl 
C 

C -  INITIALIZE  COUNTERS  - 

C 

1  =  0 
J  =  1 
PASS  =  0 
C 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  50  ITERATIONS 

C 

5  1  =  1  +  1 

IF(I.GT. J*50)  GOTO  10 
XR  =  X2-F2*(X2-X1)/(F2-F1) 

FR  =  FCN1(XR,TB,TE,HT,DEL,K,K1,K2) 

IF(PASS.EQ.O)  PRINT  100,1 
IF(PASS.EQ. 1)  PRINT  200 
IF(PASS.GE. 1)  PRINT  300,I,XR,FR 
IF(PASS.GE. 1)  PASS  =  PASS  +  1 
C 

C -  CHECK  STOPPING  CRITERIA  - 

C 

XERR  =  DABS(X1-X2)/2.0D0 

IF(XERR.LE.XTOL.OR.DABS(FR) .LE.FTOL)  RETURN 
C 

c - find  new  point - 

C 

IF(FR*F1.GE.0.0D0)  THEN 
XI  =  XR 
Fl  =  FR 

IF(FR*FSAVE.GT.0.0D0)  F2  =  F2/2.0D0 
FSAVE  =  FR 
ELSE 

X2  =  XR 
F2  =  FR 

IF(FR*FSAVE.GT.0.0D0)  Fl  =  F1/2.0D0 
FSAVE  =  FR 
ENDIF 
GOTO  5 
C 

C - FAIL  TO  FIND  ROOT,  CONTINUE  SEARCH  7 - 

C 

10  PRINT  400,J*50 
READ  500, ANS 

IF(ANS.NE. • y .AND.ANS.NE. ’N' )  THEN 
PRINT  600 
GOTO  10 
ENDIF 

IF(ANS.EQ. ’N* )  GOTO  15 
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c 

c 

c 


J  =  J  +  1 
1  =  1-1 
GOTO  5 

-  FAIL  TO  FIND  ROOT,  LOOK  AT  FUNTION  VALUES  7 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. 'Y' .AND.ANS.NE. 'N* )  THEN 
PRINT  600 
GOTO  15 
ENDIF 

IF(ANS.EQ. ‘N* )  STOP 
XI  =  XISAV 
X2  =  X2SAV 
FI  =  FISAV 
F2  =  F2SAV 
1  =  0 
J  =  1 
PASS  =  1 
GOTO  5 
C 

c  - 

c 


100  FORMAT('+',  'ITERATION  :  ',13,'  ') 

200  FORMAT (/,'  COMMENCING  SEARCH  FOR  ROOT  IN  INTERVAL....  ',/) 

300  FORMAT('  AT  ITERATION ', I3 , 3X, '  X  =  •,D15.6,3X, '  F(X)  =  ’,015.6) 

400  FORMAT(/,’  FAILED  TO  FIND  ROOT  AFTER  *,13, *  ITERATIONS.’  , 

&  ’  CONTINUE  SEARCH  (Y  OR  N)  ?’  ) 

500  FORMAT(A2) 

600  FORMAT(/’  *****************  INCORRECT  ANSWER  ********************) 
700  FORMAT(/,’  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N  )7  ’,/, 

&  •  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  ’,/, 

&  ’  IF  F(XR)  IS  OSCILLATING  ABOUT  A  CERTAIN  VALUE  OR  ITS  ’,/, 

&  ’  ABSOLUTE  VALUE  IS  INCREASING  THEN  THERE  IS  LITTLE  ’  , / , 

&  ’  PROBABILITY  THAT  A  ROOT  WILL  BE  FOUND.  ’) 


END 

C 

C3A******  RTSY:  HEADER2  ********************************************** 
C 

C  SUBROUTINE  RTSY ( TB, TE , L, HT, DEL, K,ABS, EMIS , E, QI , Q, EFF, FLAG) 

C 

C  FIN  SHAPE:  RECTANGLE 
C 

C  PROBLEM  TYPE:  SYNTHESIS 
C 

C  INPUT:  TB,L,DEL,K,ABS,EMIS,E,Q 
C 

C  OUTPUT :  HT , TE , QI , EFF 
C 

c - 

c 

C  PARAMETERS : 
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c 

C  TB  -  TEMPERATURE  AT  THE  BASE 

C  TE  -  TEMPERATURE  AT  THE  TIP 

C  HT  -  HEIGHT  OF  THE  FIN 

CL-  LENGTH  OF  THE  FIN 

C  DEL  -  THICKNESS  OF  THE  FIN 

C  K  -  CONDUCTIVITY  OF  THE  FIN 

C  ABS  -  ABSORPTIVITY  FOR  FIN 

C  EMIS  -  EMISSIVITY  OF  THE  FIN 

C  E  -  EXTERNAL  HEAT  INCIDENT  ON  THE  FIN 

C  QI  -  IDEAL  HEAT  DISSIPATED  BY  THE  FIN 

C  Q  -  REAL  HEAT  DISSIPATED  BY  THE  FIN 

C  EFF  -  EFFICIENCY  OF  THE  FIN 

C  FLAG  -  0  =  CONVERGENCE,  1  =  DIVERGENCE 

C 

C - 

c 

C  FUNCTIONS : 

C 

C  FCN2(X,TB,TE,L,DEL,K,K1,K2,Q) : 

C 

C  FUNCTION  WHOSE  ROOT  IS  FOUND  (BASE  DERIV  &  TIP  TEMP) 

C 

C - 

C 

C  SUBROUTINES ; 

C 

C  RKFSY2 (DERIV2,XO,TO,TEND,F,H,DEL,K,K1,K2,TOL) ; 

C 

C  SOLVES  SECOND  ORDER  NONLINEAR  DE  BY  RUNGE-KUTTE-FELHBERG  METHOD 
C 

C  DERIV2(T,F,DEL,K,K1,K2) : 

C 

C  COMPUTES  DERIVATIVES  FOR  RKFSY2 
C 

C  MDLIN2(FCN2,X1,X2,XR,XT0L,FT0L,NLIM,TB,TE,L,DEL,K,K1,K2,Q) : 

C 

C  FINDS  THE  ROOT  OF  FCN2  BY  THE  METHOD  OF  MODIFIED  LINEAR 
C  INTERPOLATION 
C 

C3B  *****  RTSY:  MAIN2  ************************************************ 
C 

SUBROUTINE  RTSY ( TB, TE , L, HT, DEL, K, ABS, EMIS , E , QI , Q, EFF, FLAG) 

REAL*8  TB,TE,L,HT,DEL,K,ABS,EMIS,E,QI,Q,EFF,X1,X2,XR, 

&  SB,K1,K2,FCN2,F1,F2 

INTEGER  I, J, PASS, FLAG 
CHARACTER* 2  ANS 
EXTERNAL  FCN2 


C 

C -  DEFINE  CONSTANTS 

C 

SB  =  5.67051D-12 
K1  =  2 .0D0*SB*EMIS 
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K2  =  E*ABS 
C 

C -  SCALE  INPUTS  - 

C 

DEL  =  DEL*100.0D0 

L  =  L*100.0D0 

K  =  K*0.01D0 

C 

C -  START  COUNTERS 

C 

PASS  =  0 


1  =  0 
J  =  1 
c 

C -  START  SEARCH  WITH  INTERVAL  [0.1, 0.2]  - 

C 

PRINT  200 
XI  =  O.IDO 
X2  =  0.2D0 

Fl  =  FCN2(Xl,TB,TE,L,DEL,K,Kl,K2,Q) 

5  1  =  1  +  1 
C 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  100  INTERVALS 

C 


IF(I.GT.J*100)  GOTO  10 

IF(PASS.EQ.O)  PRINT  100,1 

IF(PASS.EQ. 1)  PRINT  200 

IF(PASS.GE.l)  PRINT  300 , Xl , Fl , X2 , F2 

F2  =  FCN2(X2,TB,TE,L,DEL,K,K1,K2,Q) 

C 

C -  CHECK  FOR  FUNCTION  VALUES  OPPOSITE  IN  SIGN 

C 

IF(F1*F2.GE.0.0D0)  THEN 
XI  =  X2 
Fl  =  F2 

IF(X2.LT.1.0D0)  X2  =  X2  +  O.IDO 
IF(X2.GE.1.0D0)  X2  =  X2  +  l.ODO 
IF(PASS.GE. 1)  PASS  =  PASS  +  1 
GOTO  5 
END  IF 
C 

C - IF  INTERVAL  FOUND,  GO  FIND  ROOT - 

C 

PRINT  201 
GOTO  20 
C 

C - IF  NO  INTERVAL  FOUND,  CONTINUE  SEARCH  7  — 

C 

10  PRINT  400,J*100 
READ  500,ANS 

IF(ANS.NE. • Y' .AND.ANS.NE. ’N’ )  THEN 
PRINT  600 
GOTO  10 
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ENDIF 

IF(ANS.EQ. 'N* )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 
C 

C -  IF  UNABLE  TO  FIND  INTERVAL,  LOOK  AT  FUNCTION  VALUES  7 

C 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. 'Y' .AND.ANS.NE. 'N' )  THEN 
PRINT  600 
GOTO  15 
ENDIF 

IF(ANS.EQ. ‘N* )  THEN 
FLAG  =  1 
RETURN 
ENDIF 


PASS  =  1 
1  =  0 
J  =  1 

XI  =  O.lDO 
X2  =  0.2D0 

Fl  =  FCN2 (X1,TB,TE,L,DEL,K,K1,K2,Q) 

GOTO  5 
C 

C -  GET  ROOT  IN  INTERVAL:  FIN  HEIGHT  AND  TEMP  AT  TIP 


C 

20  CALL  MDLIN2(FCN2,X1,X2,XR,TB,TE,L,DEL,K,K1,K2,Q) 


HT  =  XR 

c 

C -  CALCULATE  IDEAL  HEAT  DISSIPATED  - 

C 

QI  =  2.0D0*SB*EMIS*HT*L*TB**4.0D0 

c 

C -  CALCULATE  FIN  EFFICIENCY  - 

C 

EFF  =  Q/QI 
C 

C -  SCALE  OUTPUTS  - 

C 

HT  =  HT*0.01D0 
DEL  =  DEL* 0.0 IDO 
L  =  L*0.01DO 
K  =  K*100.0D0 
C 

C -  FORMAT  - 

C 

100  FORMAT (•+•,' LOOKING  IN  INTERVAL:  ’,13, •  ’) 

200  FORMAT(/,’  COMMENCING  SEARCH  FOR  INTERVAL  TO  BRACKET  ROOT _  ',/) 

201  FORMAT (/,•  INTERVAL  FOUND,  WILL  NOW  LOOK  FOR  ROOT  _  ’,/) 

300  FORMATC  XL  = ' , 2X, DIO . 4 , 2X, -FCXL)  =• ,2X,D10.4,2X, 

&  ’XR  =• ,2X,D10.4,2X, 'F<XR)  =*,2X,D10.4) 


92 


400  FORMAT (/,•  NO  ROOT  FOUND  IN  ',13, ’  INTERVALS.  CONTINUE  SEARCH*  , 
&  *  (Y  OR  N)  7*  ) 

500  FORMAT(A2) 

600  FORMAT(/,'  ****************  INCORRECT  ENTRY  I  *****************  •) 
700  FORMAT (/,*  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N)7  *,/, 

&  *  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  IF  *,/, 
&  '  F(XR)  IS  OF  THE  SAME  SIGN  AS  F(XL)  AND  ITS  ABSOLUTE  VALUE 
&  *  IS  INCREASING,  THERE  IS  LITTLE  PROBABILITY  OF  FINDING  AN  *,/, 
&  *  INTERVAL.  THE  SAME  GOES  FOR  VALUES  OF  F(XR)  THAT  OSCILLATE’,/, 
&  *  AROUND  A  CERTAIN  VALUE.  FOR  F(XR)  TO  BECOME  OPPOSITE  IN  *,/, 
&  *  SIGN  TO  F(XR),  IT  MUST  PASS  THROUGH  ZERO.  *) 

RETURN 
END 
C 

C3C******  RTSY;  FCN2 (X,TB,TE,L,DEL,K,K1,K2,Q)  *********************** 

C 

C  PARAMETERS : 

C 

C  X  -  INDEPENDENT  VARIABLE  (BASE  DERIV) 

C  TB  -  BASE  TEMP 
C  TE  -  TIP  TEMP 
C  DEL  -  FIN  THICKNESS 
C  K  -  FIN  CONDUCTIVITY 
C  K1  -  CONSTANT  =  2.0*SB*EMIS 
C  K2  -  CONSTANT  =  ABS*E 
C  Q  -  HEAT  INPUT  THRU  BASE 
C 

C - 

c 

REAL*8  FUNCTION  FCN2 ( X, TB, TE,L,DEL, K, K1 , K2 , Q) 

REAL*8  X,X0,XFINAL,T0(2) ,TEND(2) ,F(2) ,H,TB,TE,DEL, 


&  K,Kl,K2,L,Q,TSAVE,TOL 

EXTERNAL  DERIV2 
C 

C -  INITIALIZE  LEFT  END  BC*S;  KNOW  FCT  VALUE  &  DERIV 


C 

T0( 1)=  TB 

T0(2)=  -Q/(K*DEL*L) 

C 

C - INITIALIZE  RIGHT  END  BC*S  TO  0 

C 

TEND(1)=0.0D0 

TEND(2)=0.0D0 

C 

C -  DEFINE  PARAMETERS  - 

C 

IF  (X.LT.1.0)  H  =  0.05D0 

IF  (X.GE.1.0)  H  =  O.IDO 

TOL  =  O.OOOIDO 

XO  =  O.ODO 

XFINAL  =  X 

TSAVE  =  1000. ODD 

10  IF(X0.LT.XFINAL+0.0001D0)  THEN 
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c 

c 

c 


CALL  RKFSY2(DERIV2,XO,TO,TEND,F,H,DEL,K,K1,K2,TOL) 


20 


C 

C3D 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c — 
c 


c 

c 

c 


-  INCREASING  T  MEANS  DIVERGENCE 

IF(TEND(1) .GE.TSAVE)  GO  TO  20 
T0(1)  =  TEND(l) 

TSAVE  =  TEND(l) 

TO (2)  =  TEND(2) 

GOTO  10 
ENDIF 
CONTINira: 

FCN2  =  TEND (2) 

TE  =  TEND(l) 

RETURN 

END 


RTSY;  RKFSY2(DERIV2,XO,TO,TEND,F,H,DEL,K,K1,K2,TOL) 
PARAMTERS : 

RKFSY2  -  SUBROUTINE  THAT  SOLVES  A  SYSTEM  OF  2  FIRST  ORDER 

DIFFERENTIAL  EQUATIONS  BY  THE  RUNGE-KUTTA-FEHLBERG 
METHOD.  THE  EQUATIONS  ARE  OF  THE  FORM: 

DT/DX  =  Y  =  Fl(X,T) 

DY/DX  =  F2(X,T,Y) 

DERIV2  -  A  SUBROUTINE  THAT  COMPUTES  VALUES  OF  THE  2  DERIVATIVES. 

XO  -  THE  INITIAL  VALUE  OF  INDEPENDENT  VARIABLE 

TO  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 

TEND  -  AN  ARRAY  THAT  RETURNS  THE  PINAL  VALUES  OF  THE  FUNCTIONS 

F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 

H  -  THE  INCREMENT  TO  T,  THE  STEP  SIZE 

DEL  -  FIN  THICKNESS 

K  -  FIN  CONDUCTIVITY 

K1  -  CONSTANT  =  2.0*SB*EMIS 

K2  -  CONSTANT  =  E*ABS 

TOL  -  TOLERANCE 

TWRK  -  AN  ARRAY  USED  TO  HOLD  INTERMEDIATE  VALUES  DURING  THE 
COMPUTATION.  IT  MUST  BE  DIMENSIONED  OF  SIZE  6X2. 


SUBROUTINE  RKFS Y2 ( DERI V2 , XO , TO , TEND , F , H , DEL , K , K 1 , K2 , TOL ) 
REAL*8  X0,T0(2) ,TEND(2) ,F(2) ,TWRK(6,2),H,K,Kl,K2,TOL,DEL, 
&  ERROR, SUM, STOREH,XEND 

INTEGER  I 

-  INITIALIZE  FOR  INTERVAL  [XO,  XO+H]  - 


XEND  =  XO+H 
STOREH  =  H 
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c 

C - CHECK  TO  SEE  IF  WE  ARE  FINISHED - 

C 

1  IF(XO.GE.XEND)  THEN 
H  =  STOREH 
RETURN 
ELSE 
C 

C - GET  FIRST  ESTIMATE  OF  THE  DELTA  X’S - 

C 

CALL  DERIV2(T0,F,DEL,K,K1,K2) 

DO  10  I  =  1,2 

TWRK(1,I)  =  H*F(I) 

TEND(I)  =  T0(1)+TWRK(1,I)/4,0D0 
10  CONTINUE 
C 

c -  get  second  estimate  - 

c 

CALL  DERIV2 (TEND, F, DEL, K,K1,K2) 

DO  20  I  =  1,2 

TWRK(2,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(TWRK(l,I)*3.ODO+TWRK(2,I)*9.OD0)/32.0D0 
20  CONTINUE 
C 

C -  GET  THIRD  ESTIMATE  - 

C 

CALL  DERIV2(TEND,F,DEL,K,K1,K2) 

DO  30  I  =  1,2 

TWRK(3,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(TWRK( 1,I)*1932.0D0-TWRK(2,I)*7200.0D0 
&  +  TWRK(3,I)*7296 .0D0)/2197 .ODO 

30  CONTINUE 
C 

C -  GET  FOURTH  ESTIMATE  - 

C 

CALL  DERIV2 (TEND, F, DEL, K,K1,K2) 

DO  40  I  =  1,2 

TWRK(4,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(439.0D0*TWRK(1,I)/216.0D0-8.0D0*TWRK(2,I) 
&  +  3680.0DO*TWRK(3,I)/513.0DO-845.0DO*TWRK(4,I)/4104.0DO) 

40  CONTINUE 
C 

c -  get  fifth  estimate  - 

c 

CALL  DERIV2( TEND, F, DEL, K,K1,K2) 

DO  50  I  =  1,2 

TWRK(5,I)  =  H*F(I) 

TEND(I)  =  TO(I)-8.0DO*TWRK(1,I)/27.0+2.0DO*TWRK(2,I) 

&  -  3544.0D0*TWRK(3,I)/2565.0D0+1859.0D0*TWRK(4,I)/4104.0D0 

&  -  11.0D0*TWRK(5,I)/40.0D0 

50  CONTINUE 
C 

c -  get  SIXTH  ESTIMATE  - 
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c 

CALL  DERIV2 (TEND, F, DEL, K,K1,K2) 

DO  60  I  =  1,2 

TWRK(6,I)  =  H*F(I) 

60  CONTINUE 

C 

C -  ESTIMATE  THE  ERROR  BY  COMPUTING  THE  DIFFERENCE  BETWEEN 

C  THE  FOURTH  AND  FIFTH  ORDER  EQUATIONS 

C 

SUM  =  O.ODO 
ERROR  =  O.ODO 
DO  70  I  =  1,2 

SUM  =  DABS(TWRK(1,I)/360.0D0-128.0D0*TWRK(3,I)/4275.0D0 

&  -  2197.0D0*TWRK(4,I)/75240.0D0+TWRK(5,I)/50.0D0 

&  +  2.0D0*TWRK(6,I)/55.0D0) 

ERROR  =  DMAXl (ERROR, SUM) 

70  CONTINUE 

C 

c -  IF  ERROR  IS  LESS  THAN  TOLERANCE,  COMPUTE  X  AT 

C  END  OF  THE  INTERVAL  FROM  A  WEIGHTED  AVERAGE  OF  THE  SIX 

C  ESTIMATES  &  THEN  RETURN 

C 

IF ( ERROR. LT.TOL)  THEN 
DO  80  I  =  1,2 


TEND(I)  =  T0(I)+16.0D0*TWRK(1,I)/135.0D0 
&  +  6656.0d0*TWRK(3,I)/12825.0D0 

&  +  28561. 0D0*TWRK(4,I)/56430. 0-9. 0D0*TWRK(5,I)/50.0D0 

&  +  2.0D0»TWRK(6,I)/55.CD0 

T0(I)  =  TEND(I) 

80  CONTINUE 

XO  =  XO+H 
END  IF 
C 

C -  IF  ERROR  IS  GREATER  THAN  TOLERANCE,  DECREASE  STEP  &  - 

C  GO  AGAIN 

C 

IF ( ERROR. GT.TOL)  THEN 
H  =  H/2.0D0 
ENDIF 
C 

C -  IF  ERROR  IS  SIGNIFICANTLY  LESS  THAN  TOLERANCE,  RELAX  STEP  — 

C 

IF ( ERROR. LT.H*TOL/l O.ODO)  THEN 
H  =  H*2.0D0 
ENDIF 
C 

C -  IF  OVERSHOOT,  REDUCE  STEP  - 


C 

IF(X0  +  H.GT.XEND)  THEN 
H  =  XEND  -  XO 
ENDIF 
C 

ENDIF 
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c 

c 

GO  TO  1 
END 
C 

C3E******  RTSY;  DERIV2 (T, F,DEL, K, Kl,K2 )  ***************************** 
C 

C  PARAMETERS : 

C 

C  T  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 

C  F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 

C  DEL  -  FIN  THICKNESS 

C  K  -  FIN  CONDUCTIVITY 

C  K1  -  CONSTANT  =  2 . 0*SB*EMIS 

C  K2  -  CONSTANT  =  E*ABS 

C 

C - 

c 

SUBROUTINE  DERI V2 ( T , F , DEL , K , K1 , K2 ) 

REALMS  T(2) ,F(2) ,K,K1,K2,DEL 
F(l)  =  T(2) 

F(2)  =  K1*T( 1)»*4.0D0/(K»DEL)  -  K2/(K*DEL) 

RETURN 

END 

C 

C3F  *****  RTSY;  MDLIN2 (FCN2 , XI , X2 , XR, XTOL, FTOL,NLIM, TB, TE,L, DEL, 

C  K,K1,K2,Q) 

C 

C  PARAMETERS : 

C 

C  FCN2  -  FUNCTION  THAT  COMPUTES  VALUES  FOR  F.  MUST  BE  DECLARED 
C  EXTERNAL  IN  CALLING  PROGRAM 

C  Xl,X2  -  INITIAL  VALUES  OF  X.  F(X)  MUST  CHANGE  SIGNS  AT  THESE 
C  POINTS 

C  XR  -  RETURNS  THE  ROOT  TO  THE  MAIN  PROGRAM 

C  XTOL  -  TOLERANCE  FOR  X 

C  FTOL  -  TOLERANCE  FOR  F 

C  NLIM  -  LIMIT  TO  NUMBER  OF  ITERATIONS 

C  TB  -  BASE  TEMPERATURE 

C  TE  -  TIP  TEMPERATURE 

CL-  FIN  LENGTH 

C  DEL  -  FIN  THICKNESS 

C  K  -  FIN  CONDUCTIVITY 

C  K1  -  CONSTANT  =  2.0*SB*EMIS 

C  K2  -  CONSTANT  =  E*ABS 

C  Q  -  HEAT  INPUT  THRU  BASE 

C 

C - 

C 

SUBROUTINE  MDLIN2 ( FCN2 , X 1 , X2 , XR, TB , TE , L , DEL , K , Kl , K2 , Q ) 

REAL*8  XERR,FSAVE,Fl,F2,FR,Xl,X2,XR,XTOL,FTOL,TB,TE,L,DEL, 

&  K, Kl , K2 , Q, FCN2 , XISAV, X2SAV, FlSAV, F2SAV 

INTEGER  I, J, PASS 
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CHARACTER* 2  ANS 

DATA  XTOL,FTOL/0. 000 100,0.0000100/ 

C 

c -  INITIALIZE  VALUES  - 

C 

PRINT  200 
XlSAV  =  XI 
X2SAV  =  X2 

Fl  =  FCN2 (X1,TB,TE,L,DEL,K,K1,K2,Q) 

F2  =  FCN2(X2,TB,TE,L,DEL,K,K1,K2,Q) 

FlSAV  =  Fl 
F2SAV  =  F2 
FSAVE  =  Fl 
C 

C -  INITIALIZE  COUNTER  - 

C 

1  =  0 
J  =  1 
PASS  =  0 
c 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  50  ITERATIONS 

C 

5  1  =  1  +  1 

IF(I.GT.J*50)  GOIO  10 
XR  =  X2-F2*(X2-X1)/(F2-F1) 

FR  =  FCN2(XR,TB,TE,L,DEL,K,K1,K2,Q) 
IF(PASS.EQ.O)  PRINT  100,1 
IF(PASS.EQ.l)  PRINT  200 
IF(PASS.GE. 1)  PRINT  300,I,XR,FR 
IF(PASS.GE.l)  PASS  =  PASS  +  1 
C 

C -  CHECK  STOPPING  CRITERIA  - 

C 

XERR  =  DABS(X1-X2)/2.0D0 

IF(XERR.LE.XTOL.OR.DABS(FR) .LE.FTOL)  RETURN 
C 

c - find  new  point - 

c 

IF(FR*F1,GE. 0.000)  THEN 
XI  =  XR 
Fl  =  FR 

IF(FR»FSAVE.GT.C.0D0)  F2  =  F2/2.0D0 
FSAVE  =  FR 
ELSE 
X2  =  XR 
F2  =  FR 

IF (FR*FSAVE.GT. 0.000)  Fl  =  F1/2.0D0 
FSAVE  =  FR 
ENDIF 
GOTO  5 
C 

C - FAIL  TO  FIND  ROOT,  CONTINUE  SEARCH  ? - 

C 
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c 

c- 

c 


10  PRINT  400,J*50 
READ  500, ANS 

IF(ANS.NE. ’Y* .AND.ANS.NE. ’N* )  THEN 
PRINT  600 
GOTO  10 
END  IF 

IF(ANS.EQ. 'N* )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 


FAIL  TO  FIND  ROOT,  LOOK  AT  FUNTION  VALUES  ? 


15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. ’ y .AND.ANS.NE. 'N' )  THEN 
PRINT  600 
GOTO  15 


ENDIF 

IF(ANS.EQ.  ' 
XI  =  XISAV 
X2  =  X2SAV 
FI  =  FlSAV 
F2  =  F2SAV 
1  =  0 
J  =  1 
PASS  =  1 
GOTO  5 


N • )  STOP 


C 

C 


100  FORMAT('+',  'ITERATION  :  ',13,'  •) 

200  FORMAT(/,'  COMMENCING  SEARCH  FOR  ROOT  IN  INTERVAL _  ',/) 

300  FORMAT('  AT  ITERATION ', 13 , 3X, '  X  =  •,D15.6,3X,'  F(X)  =  ',015.6) 

400  FORMAT(/,'  FAILED  TO  FIND  ROOT  AFTER  ’,13,'  ITERATIONS.’  , 

&  '  CONTINUE  SEARCH  (Y  OR  N)  7’  ) 

500  FORMAT (A2) 

600  FORMAT(/'  *****************  INCORRECT  ANSWER  *******************’) 
700  FORMAT (/,'  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N  )?  ',/, 

&  '  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  ’,/, 

&  '  IF  F(XR)  IS  OSCILLATING  ABOUT  A  CERTAIN  VALUE  OR  ITS  ’,/, 

&  ’  ABSOLUTE  VALUE  IS  INCREASING  THEN  THERE  IS  LITTLE  ',/, 

&  '  PROBABILITY  THAT  A  ROOT  WILL  BE  FOUND.  ’) 

END 


C4A  ******  TPAN;  HEADER3*’************»**»***»»»**»******************»* 
C 

C  SUBROUTINE  TPAN ( TB , TE , L , HT , DEL 0 , DELE ,K,ABS,EMIS,E,QI,Q, EFF , FLAG ) 

C 

C  FIN  SHAPE:  TRAPEZOID 
C 

C  PROBLEM  TYPE:  ANALYSIS 
C 
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C  INPUT;  TB,L,HT,DELO,DELE,K,ABS,EMIS,E 
C 

C  OUTPUT;  Q,TE,QI,EFF 
C 

C - 

c 

C  PARAMETERS ; 

C 

C  TB  -  TEMPERATURE  AT  THE  BASE  OF  THE  FIN 

C  TE  -  TEMPERATURE  AT  THE  TIP  OF  THE  FIN 

C  HT  -  HEIGHT  OF  THE  FIN 

CL-  LENGTH  OF  THE  FIN 

C  DELO  -  THICKNESS  AT  FIN  BASE 

C  DELE  -  THICKNESS  AT  FIN  TIP 

C  K  -  CONDUCTIVITY  OF  THE  FIN 

C  ABS  -  ABSORPTIVITY  OF  THE  FIN 

C  EMIS  -  EMISSIVITY  OF  THE  FIN 

C  E  -  EXTERNAL  HEAT  INCIDENT  ON  THE  FIN 

C  QI  -  IDEAL  HEAT  DISSIPATED  BY  THE  FIN 

C  Q  -  REAL  HEAT  DISSIPATED  BY  THE  FIN 

C  EFF  -  EFFICIENCY  OF  THE  FIN 

C  FLAG  -  0  =  CONVERGENCE,  1  =  DIVERGENCE 

C 

C - 

c 

C  FUNCTIONS ; 

C 

C  FCN3(X, TB,TE,HT, DELO, DELE, K, Kl, K2) ; 

C 

C  FUNCTION  WHOSE  ROOT  IS  FOUND  (BASE  DERIV  &  TIP  TEMP) 

C 

C - 

c 

C  SUBROUTINES; 

C 

C  RKFSY3(DERIV3,X0, TO, TEND, F,H,HT, DELO, DELE, K,Kl,K2,TOL) ; 

C 

C  SOLVES  A  SECOND  ORDER  NONLINEAR  DE  BY  RUNGE-KUTTE-FEHLBERG  METHOD 
C 

C  DERIV3(X,T,F,HT, DELO, DELE, K,K1,K2) ; 

C 

C  COMPUTES  DERIVATIVES  FOR  RKFSY3 
C 

C  MDLIN3 ( FCN3 , X 1 , X2 , XR , XTOL , FTOL , NLIM , TB , TE , HT , DEL  0 , DELE , K , K 1 , K2 ) ; 

C 

C  FINDS  THE  ROOT  OF  EQUATION  FCN3  BY  THE  METHOD  OF  MODIFIED 
C  LINEAR  INTERPOLATION 
C 

C4B  *****  TPAN;  MAIN3***********»******»»***************************** 

C 
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SUBROUTINE  TPAN (TB, TE, L, HT, DELO , DELE, K,ABS , EMIS , E, QI , Q, EFF, FLAG) 
REAL*8  TB,TE,L,HT,DEL0,DELE,K,ABS,EMIS,E,QI,Q,EFF,X1,X2,XR, 

&  SB,TBDER,K1,K2,F1,F2,FCN3 

INTEGER  I, J, PASS, FLAG 
CHARACTER* 2  ANS 
EXTERNAL  FCN3 
C 

C - DEFINE  CONSTANTS - - - 

C 

SB  =  5.67051D-12 
Kl  =  2.0D0*SB*EMIS 
K2  =  E*ABS 
C 

C -  SCALE  INPUTS  - 

C 

HT  =  HT*100.0D0 
DELO  =  DEL0*100.0D0 
DELE  =  DELE*100.0D0 
L  =  L*100.0D0 
K  =  K*0.01D0 
C 

C -  START  COUNTERS  - 

C 

PASS  =  0 


1  =  0 
J  =  1 
c 

c -  START  SEARCH  WITH  INTERVAL  [-1,0]  - 

C 

PRINT  200 
XI  =  O.ODO 
X2  =  -l.ODO 

FI  =  FCN3(X1,TB,TE,HT, DELO, DELE, K,K1,K2) 

5  1  =  1  +  1 
C 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  100  INTERVALS 

C 


IF(I.GT.J»100)  GOTO  10 
IF(PASS.EQ.O)  PRINT  100,1 
IF(PASS.EQ. 1)  PRINT  200 
IF(PASS.GE, 1)  PRINT  300 , Xl , Fl , X2 , F2 


F2  =  FCN3(X2,TB,TE,HT, DELO, DELE, K,K1,K2) 

C 

C -  CHECK  FOR  FUNCTION  VALUES  OPPOSITE  IN  SIGN 

C 


IF(F1*F2.GE. O.ODO)  THEN 
XI  =  X2 
Fl  =  F2 

X2  =  X2  -  l.ODO 
IF(PASS.GE.l)  PASS  =  PASS  +  1 
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r>  n 


GOTO  5 
ENDIF 
C 

C - IF  INTERVAL  FOUND,  GO  FIND  ROOT - 

C 

PRINT  201 
GOTO  20 
C 

C -  IF  NO  INTERVAL  FOUND,  CONTINUE  SEARCH  7  - 

C 

10  PRINT  400,J»100 
READ  500, ANS 

IF(ANS.NE. 'y .AND.ANS.NE. 'N* )  THEN 
PRINT  600 
GOTO  10 
ENDIF 

IF(ANS.EQ. 'N* )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 
C 

C - IF  UNABLE  TO  FIND  INTERVAL,  LOOK  AT  FUNCTION  VALUES  ? 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. 'y .AND.ANS.NE. ’N* )  THEN 
PRINT  600 
GOTO  15 
ENDIF 

IF(ANS.EQ. ’N’ )  THEN 
FLAG  =  1 
RETURN 
ENDIF 
PASS  =  1 
1  =  0 
J  =  1 

XI  =  O.ODO 
X2  =  -l.ODO 

FI  =  FCN3(X1,TB,TE,HT,DEL0,DELE,K,K1,K2) 

GOTO  5 

-  GET  ROOT  IN  INTERVAL:  DERIV  AT  BASE  &  TEMP  AT  TIP  — 

20  CALL  MDLIN3(FCN3,X1,X2,XR,TB,TE,HT,DEL0,DELE,K,K1,K2) 
TBDER  =  XR 


C 

C -  CALCULATE  IDEAL  HEAT  DISSIPATED 

C 

01  =  2.0D0*SB*EMIS*HT*L*TB**4.0D0 
C 
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c - CALCUI*ATE  REAL*  8  HEAT  DISSIPATED - 

C 

Q  =  -K*DELO*L*TBDER 
C 

C -  CALCULATE  FIN  EFFICIENCY  - 

C 

EFF  =  Q/QI 
C 

C -  SCALE  OUTPUTS  - - 

C 

HT  =  HT*0.01D0 
DELO  =  DEL0*0.01D0 
DELE  =  DELE* 0.0 IDO 
L  =  L*0.01D0 
K  =  K*100.0D0 
RETURN 
C 

C -  FORMAT  - 

C 

100  FORMAT (•+•,• LOOKING  IN  INTERVAL;  *,13, *  ') 

200  FORMAT(/,'  COMMENCING  SEARCH  FOR  INTERVAL  TO  BRACKET  ROOT _  ',/) 

201  FORMAT (/,'  INTERVAL  FOUND,  WILL  NOW  LOOK  FOR  ROOT  ....  *,/) 

300  FORMAT('  XL  = • , 2X, DlO . 4 , 2X, • F(XL)  = • , 2X, DIO . 4, 2X, 

&  'XR  =• ,2X,D10.4,2X, •F(XR)  =',2X,D10.4) 

400  FORMAT (/,'  NO  ROOT  FOUND  IN  *,13, •  INTERVALS.  CONTINUE  SEARCH’  , 
&  •  (Y  OR  N)  ?•  ) 

500  FORMAT (A2) 

600  FORMAT{/,’  ****************  INCORRECT  ENTRY  I  *****************  •) 
700  FORMAT (/,’  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N)? 

&  •  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  IF 
&  •  F(XR)  IS  OF  THE  SAME  SIGN  AS  F{XL)  AND  ITS  ABSOLUTE  VALUE 
&  '  IS  INCREASING,  THERE  IS  LITTLE  PROBABILITY  OF  FINDING  AN  ’,/, 
&  •  INTERVAL.  THE  SAME  GOES  FOR  VALUES  OF  F(XR)  THAT  OSCILLATE',/, 
&  •  AROUND  A  CERTAIN  VALUE.  FOR  F(XR)  TO  BECOME  OPPOSITE  IN  ’,/, 
&  •  SIGN  TO  F(XR),  IT  MUST  PASS  THROUGH  ZERO.  ') 

END 
C 

C4C******  TPAN:  FCN3 ( X,TB, TE, HT, DELO, DELE, K, K1 ,K2 )  ****************** 

C 

C  PARAMETERS : 

C 

C  X  -  INDEPENDENT  VARIABLE  (BASE  DERIV) 

C  TB  -  BASE  TEMPERATURE 

C  TE  -  TIP  TEMPERATURE 

C  HT  -  FIN  HEIGHT 

C  DELO  -  FIN  THICKNESS  AT  BASE 

C  DELE  -  FIN  THICKNESS  AT  TIP 

C  K  -  FIN  CONDUCTIVITY 
C  Kl  -  CONSTANT  =  2.0*SB*EMIS 
C  K2  -  CONSTANT  =  ABS*E 
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c 

c - 

c 

REAL»8  FUNCTION  FCN3 (X, TB, TE, HT, DELO, DELE, K, Kl , K2 ) 

REALMS  X,X0,XFINAL,T0(2) ,TEND(2) ,F(2) ,H,TB,TE,HT,DELO,DELE, 
&  K,K1,K2,TSAVE,T0L 

EXTERNAL  DERIV3 
C 

C -  INITIALIZE  LEFT  END  POINT  BC*S;  KNOW  FCT  VALUE,  GUESS 

C  DERIVATIVE 


C 

T0(1)=  TB 
T0(2)=  X 
C 

C - INITIALIZE  RIGHT  END  POINT  BC*S  TO  0 - 

C 

TEND( 1)=0.0D0 
TEND(2)=O.ODO 
C 

C -  DEFINE  PARAMETERS  - 

C 

H=  O.IDO 
TOL  =  O.OOOIDO 
XO  =  O.ODO 
XFINAL  =  HT 
TSAVE  =  1000. ODO 

10  IF(X0.LT.XFINAL+0.0001D0)  THEN 

CALL  RKFSy3 ( DERI V3 , XO , TO , TEND , F, H, HT, DELO , DELE , X, Kl , K2 , TOL ) 

C 

C -  INCREASING  T  MEANS  DIVERGENCE  - 

C 

IF  (TEND( 1) .GE. TSAVE)  GOTO  20 
T0(1)  =  TEND(l) 

TSAVE  =  TEND(l) 

TO (2)  =  TEND(2) 

GO  TO  10 
END  IF 

20  CONTINUE 

FCN3  =  TEND(2) 

TE  =  TEND(l) 

RETURN 

END 

C 

C4D  *****  TPAN:  RKFSy3<DERIV3,X0, TO, TEND, P,H,HT, DELO, DELE, K,K1,K2, TOL) 
C 

C  PARAMTERS : 

C 

C  RKFSy3  -  SUBROUTINE  THAT  SOLVES  A  SySTEM  OF  2  FIRST  ORDER 
C  DIFFERENTIAL  EQUATIONS  BY  THE  RUNGE-KUTTA-FEHLBERG 

C  METHOD.  THE  EQUATIONS  ARE  OF  THE  FORM; 
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c 

C  DT/DX  =  Y  =  Fl(X,T) 

C  DY/DX  =  F2(X,T,Y) 

C 

C  DERIV3  -  A  SUBROUTINE  THAT  COMPUTES  VALUES  OF  THE  2  DERIVATIVES. 

C 

C  XO  -  THE  INITIAL  VALUE  OF  INDEPENDENT  VARIABLE 
C  H  -  THE  INCREMENT  TO  T,  THE  STEP  SIZE 

C  TO  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 
C  TEND  -  AN  ARRAY  THAT  RETURNS  THE  FINAL  VALUES  OF  THE  FUNCTIONS 
C  H  -  STEP  SIZE 
C  HT  -  HEIGHT  OF  FIN 

C  F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 
C  DELO  -  THICKNESS  AT  BASE 
C  DELE  -  THICKNESS  AT  TIP 
C  K  -  FIN  CONDUCTIVITY 
C  K1  -  CONSTANT  =  2.0*SB*EMIS 
C  K2  -  CONSTANT  =  E*ABS 
C  TOL  -  TOLERANCE 

C  TWRK  -  AN  ARRAY  USED  TO  HOLD  INTERMEDIATE  VALUES  DURING  THE 
C  COMPUTATION.  IT  MUST  BE  DIMENSIONED  OF  SIZE  6X2. 

C 

- - 

c 

SUBROUTINE  RKFS Y 3 ( DERI V3 , X 0 , TO , TEND , F , H , HT , DEL 0 , DELE , K , K 1 , K2 , TOL ) 
REAL* 8  X0,T0(2) ,TEND(2) ,F(2) ,TWRK<6,2) ,H,HT,DEL0,DELE,K,K1,K2, 


&  TOL, ERROR, SUM, STOREH,XEND 

INTEGER  I 
C 

c - INITIALIZE  FOR  INTERVAL  [XO,  XO+H] 

C 


XEND  =  XO+H 
STOREH  =  H 
C 


C - CHECK  TO  SEE  IF  WE  ARE  FINISHED 

C 

1  IF(XO.GE.XEND)  THEN 


H  =  STOREH 
RETURN 
ELSE 
C 

C - get  first  ESTIMATE  OF  THE  DELTA  X’S  - 

C 

CALL  DERIV3 { XO , TO , F, HT, DELO , DELE, K, K1 , K2 ) 
DO  10  I  =  1,2 

TWRK( 1,1)  =  H*F(I) 

TEND(I)  =  T0(I)+TWRK(1,I)/4.0D0 
10  CONTINUE 
C 

c -  get  SECOND  ESTIMATE  - 
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so  fj 


c 


CALL  DERIV3 ( XO+H/4 . ODO , TEND , F, BT, DELO , DELE , K, K1 , K2 ) 

DO  20  I  =  1,2 

TWRK(2,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(TWRK(1,I)*3.0D0+TWRK{2,I)*9.0D0)/32.0D0 
20  CONTINUE 
C 

c -  GET  THIRD  ESTIMATE  - 

C 

CALL  DERIV3 ( XO+3 , 0D0*H/ 8 . ODO , TEND, F, HT, DELO , DELE, K, Kl , K2 ) 

DO  30  I  =  1,2 

TWRK(3,I)  =  H*F(I) 

TEND(I)  =  TO(I)+(TWRK(1,I)*1932.0DO-TWRK(2,I)*7200.0DO 
4  +  TWRK(3,I)*7296.0D0)/2197.0D0 

30  CONTINUE 
C 

c - GET  FOURTH  ESTIMATE - 

C 

CALL  DERIV3(X0+12 .0*H/ 13.0, TEND, F,HT, DELO, DELE, K,K1,K2) 

DO  40  I  =  1,2 

TWRK(4,I)  =  H*F(I) 

TEND(l)  =  T0(I)+(439.0D0*TWRK(1,I)/216.0D0-8.0D0*TWRK(2,I) 
4  +  3680. 0D0*TWRK(3,I)/513.0D0-845.0D0*TWRK(4,I)/4104. ODO) 

40  CONTINUE 
C 

c -  GET  FIFTH  ESTIMATE  - 

C 

CALL  DERIV3 ( XO+H , TEND , F , HT , DELO , DELE , K, Kl , K2 ) 

DO  50  I  *  1,2 

TWRK(5,I)  =  H*F(I) 

TEND(l)  =  T0(I)-8.0D0*TWRK(1,I)/27.0+2.0D0*TWRK(2,I) 

4  -  3544. 0D0*TWRK(3,I)/2565.0D0+1859.0D0*TWRK(4,I)/4104. ODO 

4  -  11.0D0*TWRK(5,I)/40.0D0 

50  CONTINUE 

C 

C -  GET  SIXTH  ESTIMATE  - 

C 

CALL  DERIV3 ( XO+H/2 . ODO , TEND , F , HT, DELO , DELE , K, Kl , K2 ) 

DO  60  I  =  1,2 

TWRK(6,I)  =  H*F(I) 

0  CONTINUE 


-  ESTIMATE  ERROR  BY  COMPUTING  DIFFERENCE  BETWEEN  - 

THE  FOURTH  AND  FIFTH  ORDER  EQUATIONS 

SUM  =  O.ODO 
ERROR  =  O.ODO 
DO  70  I  =  1,2 

SUM  =  DABS(TWRK(1,1)/360.0D0-128.0D0*TWRK(3,I)/4275.0D0 
-  2197. 0D0*TWRK(4, 1) /75240.0D0+TWRK<5, I) /50. ODO 
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& 


+  2.0D0*TWRK(6,I)/55.0D0) 

ERROR  »  OMAXl( ERROR, SUM) 

7  0  CONTINUE 

C 

C -  IF  ERROR  LESS  THAN  TOLERANCE,  COMPUTE  X  AT  THE  END  OF 

C  INTERVAL  FROM  A  WEIGHTED  AVERAGE  OF  THE  SIX  ESTIMATES  & 

C  THEN  RETURN 

C 


IF ( ERROR. LT.TOL)  THEN 
DO  80  I  =  1,2 

TEND(I}  =  T0(I)+16.0D0*TWRK(1,I)/135.0D0 
&  +  6656.0d0*TWRK(3,I)/12825.OD0 

&  +  28561. 0D0*TWRK(4,I)/56430. 0-9. 0D0*TWRK(5,I)/50.0D0 

&  +  2.0D0*TWRK(6,1)/55.0D0 

T0(I)  =  TEND(I) 

80  CONTINUE 
XO  =  XO+H 


END  IF 
C 

C -  IF  ERROR  GREATER  THAN  TOLERANCE,  THEN  REDUCE  STEP  - 

C 

IF ( ERROR. GT.TOL)  THEN 
H  =  H/2.0D0 
END  IF 
C 

c -  IF  ERROR  IS  SIGNIFICANTLY  LESS  THAN  TOLERANCE,  RELAX  STEP 

C 

IF ( ERROR. LT.H*TOL/10.0D0)  THEN 
H  =  H*2.0D0 
ENDIF 
C 

C -  IF  OVERSHOOT,  REDUCE  STEP  - 

C 

IF(X0+H.GT.XEND)  THEN 
H  =  XEND-XO 
ENDIF 
C 

ENDIF 

C 


GO  TO  1 
END 
C 

C4E******  TPAN:  DERIV3 (X, T, F, HT, DELO, DELE,K, K1 ,K2 )  ****************** 
C 

C  PARAMETERS : 

C 

C  X  -  INDEPENDENT  VARIABLE 

C  T  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 
C  F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 
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C  DELO  -  FIN  THICKNESS  AT  BASE 
C  DELE  -  FIN  THICKNESS  AT  TIP 
C  K  -  FIN  CONDUCTIVITY 
C  Kl  -  CONSTANT  =  2.0*SB*EMIS 
C  K2  -  CONSTANT  =  E»ABS 
C 

C - 

c 

SUBROUTINE  DERI V3 ( X , T , F , HT , DELO , DELE , K , K 1 , K2 ) 

REAL*8  X,T(2) ,F(2) ,HT, DELO, dele, K,K1,K2, del, DELP 
DEL{X)  =  DEL0+( DELE -DELO )*X/HT 
DELP  *  (DELE-DELO)/HT 
F(l)  =  T(2) 

F(2)  =  -(DELP*T{2) )/DEL(X)+(K1*T(1)**4.0D0-K2)/(K*DEL(X) ) 
RETURN 
END 
C 

C4F  *****  TPAN:MDLIN3(FCN3,X1,X2,XR,XTOL,FTOL,NLIM,TB,TE,HT,DELO, 

C  DELE,K,K1,K2) 

C 

C  PARAMETERS : 

C 

C  FCN3  -  FUNCTION  THAT  COMPUTES  VALUES  FOR  F.  MUST  BE  DECLARED 
C  EXTERNAL  IN  CALLING  PROGRAM 

C  Xl,X2  -  INITIAL  VALUES  OF  X.  F(X)  MUST  CHANGE  SIGNS  AT  THESE 
C  POINTS 

C  XR  -  RETURNS  THE  ROOT  TO  THE  MAIN  PROGRAM 
C  XTOL  -  TOLERANCE  FOR  X 
C  FTOL  -  TOLERANCE  FOR  F 
C  NLIM  -  LIMIT  TO  NUMBER  OF  ITERATIONS 
C  TB  -  BASE  TEMPERATURE 
C  TE  -  TIP  TEMPERATURE 
C  HT  -  FIN  HEIGHT 
C  DELO  -  BASE  THICKNESS 
C  DELE  -  TIP  THICKNESS 
C  K  -  FIN  CONDUCTIVITY 
C  Kl  -  CONSTANT  =  2.0*SB*EMIS 
C  K2  -  CONSTANT  =  E*ABS 
C 

C - 

c 

SUBROUTINE  MDLIN3 ( FCN3 , Xl , X2 , XR, TB , TE , HT , DELO , DELE , K , Kl , K2 ) 
REAL*8  XERR, FSAVE, Fl , F2 , FR, XI , X2 , XR, XTOL, FTOL, TB, TE, HT, DELO , 
&  DELE , K, Kl , K2 , FCN3 , XlSAV, X2SAV, FISAV, F2SAV 

INTEGER  I, J, PASS 
CHARACTER* 2  ANS 

DATA  XTOL, FTOL/0. 000 IDO, O.OOOOIDO/ 

C 

C -  INITIALIZE  VALUES  - 

C 
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PRINT  200 
XlSAV  =  XI 
X2SAV  =  X2 

Fl  =  FCN3(Xl,TB,TE,HT,DEL0,DELE,K,Kl,K2) 

F2  =  FCN3(X2,TB,TE,HT,DEL0,DELE,K,K1,K2) 

FISAV  =  Fl 
F2SAV  =  F2 
FSAVE  =  Fl 
C 

C -  INITIALIZE  COUNTER  - 

C 

1  =  0 
J  =  1 
PASS  =  0 
C 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  50  ITERATIONS 

C 

5  1  =  1  +  1 

IF(I.GT.J*50)  GOTO  10 
XR  =  X2-F2*(X2-X1)/(F2-F1) 

FR  =  FCN3(XR,TB,TE,HT,DEL0,DELE,K,K1,K2) 
IF<PASS.EQ.O)  PRINT  100,1 
IF(PASS.EQ. 1)  PRINT  200 
IF(PASS.GE.l)  PRINT  300,I,XR,FR 
IF(PASS.GE.l)  PASS  =  PASS  +  1 
C 

C -  CHECK  STOPPING  CRITERIA  - 

C 

XERR  =  DABS(X1-X2)/2.0D0 

IF(XERR.LE.XTOL.OR.DABS(FR) .LE.FTOL)  RETURN 
C 

C - FIND  NEW  POINT - 

C 

IF(FR*F1.GE.0.0D0)  THEN 
XI  =  XR 
Fl  =  FR 

IF  (FR*FSAVE.GT.0.0D0)  F2  =  F2/2.0D0 
FSAVE  =  FR 
ELSE 
X2  =  XR 
F2  =  FR 

IF(FR*FSAVE.GT.0.0D0)  Fl  =  F1/2.0D0 
FSAVE  =  FR 
END  IF 
GOTO  5 
C 

C - FAIL  TO  FIND  ROOT,  CONTINUE  SEARCH  ? - 

C 

10  PRINT  400,J*50 
READ  500, ANS 
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IF(ANS.NE. • Y' .AND.ANS.NE. ’N* )  THEN 
PRINT  600 
GOTO  10 
ENDIF 

IF(ANS.EQ. 'N' )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 
C 

c -  fail  to  find  root,  look  at  function  values  ? 

c 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. 'Y* .AND.ANS.NE. *N' )  THEN 
PRINT  600 
GOTO  15 
ENDIF 

IF(ANS.EQ. ’N' )  STOP 
XI  =  XlSAV 
X2  =  X2SAV 
FI  =  FISAV 
F2  =  F2SAV 
1  =  0 
J  =  1 
PASS  =  1 
GOTO  5 
C 

C  - 

c 


100  FORMAT(*+’,  'ITERATION  :  ',13,'  •) 

200  FORMAT (/,'  COMMENCING  SEARCH  FOR  ROOT  IN  INTERVAL....  ',/) 

300  FORMATC  AT  ITERATION ', 13 , 3X, ’  X  =  •,D15.6,3X,'  F(X)  =  *,015.6) 

400  FORMAT (/,'  FAILED  TO  FIND  ROOT  AFTER  ',13,'  ITERATIONS.'  , 

&  '  CONTINUE  SEARCH  (Y  OR  N)  ?’  ) 

500  FORMAT (A2) 

600  FORMAT(/'  *****************  INCORRECT  ANSWER  *******************') 
700  FORMAT (/,'  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N  )?  ',/, 

&  '  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  ',/, 

&  '  IF  F(XR)  IS  OSCILLATING  ABOUT  A  CERTAIN  VALUE  OR  ITS  ',/, 

&  '  ABSOLUTE  VALUE  IS  INCREASING  THEN  THERE  IS  LITTLE  ',/, 

&  '  PROBABILITY  THAT  A  ROOT  WILL  BE  FOUND .  '  ) 


END 

C 

C5A  *****  TPSY  HEADER4  *********************************************** 
C 

C  SUBROUTINE  TPSY ( TB, TE , L, HT, DELO , DELE, K, ABS , EMIS , E , QI , Q, EFF, FLAG) 

C 

C  FIN  SHAPE;  TRAPEZOID 
C 

C  PROBLEM  TYPE:  SYNTHESIS 
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c 

C  INPUT;  TB,L,DELO,DELE,K,ABS,EMIS,E,Q 
C 

C  OUTPUT:  TE,HT,QI,EFF 
C 

C - 

c 

C  PARAMETERS ; 

C 

C  TB  -  TEMPERATURE  AT  THE  BASE  OF  THE  FIN 

C  TE  -  TEMPERATURE  AT  THE  TIP  OF  THE  FIN 

CL-  LENGTH  OF  THE  FIN 
C  HT  -  HEIGHT  OF  THE  FIN 
C  DELO  -  THICKNESS  AT  FIN  BASE 

C  DELE  -  THICKNESS  AT  FIN  TIP 

C  K  -  CONDUCTIVITY  OF  THE  FIN 
C  ABS  -  ABSORPTIVITY  OF  THE  FIN 
C  EMIS  -  EMISSIVITY  OF  THE  FIN 
C  E  -  EXTERNAL  HEAT  INCIDENT  ON  THE  FIN 
C  QI  -  IDEAL  HEAT  DISSIPATED  BY  THE  FIN 
C  Q  -  REAL  HEAT  DISSIPATED  BY  THE  FIN 
C  EFF  -  EFFICIENCY  OF  THE  FIN 
C  FLAG  -  0  =  CONVERGENCE,  1  =  DIVERGENCE 
C 

C - 

c 

C  FUNCTIONS:  FCN4(X,TB,TE,L,HT, DELO, DELE, K,Kl,K2,Q) 

C 

C  FUNCTION  WHOSE  ROOT  IS  FOUND  (BASE  DERIV  &  TIP  TEMP) 

C 

C - 

c 

C  SUBROUTINES : 

C 

C  MDLIN4 ( FCN4 , Xl , X2 , XR, XTOL, FTOL, NLIM, TB, HT, DEL, K, K1 , K2 ) ; 

C 

C  FINDS  THE  ROOT  OF  AN  EQUATION  FCN4  BY  THE  METHOD  OF  MODIFIED 
C  LINEAR  INTERPOLATION 
C 

C  RKFSY4{DERIV4,XBEGIN,H, TO, TEND, F,HT, DELO, DELE, K,Kl,K2,TOL) : 

C 

C  SOLVES  SECOND  ORDER  NONLINEAR  DE  BY  RUNGE-KUTTE-FEHLBERG  METHOD 
C 

C  DERIV4 (X,T,F,HT, DELO, DELE, K,K1,K2) : 

C 

C  COMPUTES  DERIVATIVES  FOR  RKFSY4 
C 

C5B  *****  TPSY;  MAIN4  ************************************************ 
C 

SUBROUTINE  TPSY ( TB, TE , L, HT , DELO , DELE, K, ABS , EMIS , E , QI , Q, EFF, FLAG) 


111 


REAL* 8  TB,TE,L,HT,DEL0,DELE,K,ABS,EMIS,E,QI,Q,EFF,X1,X2,XR, 
&  SB,K1,K2,F1,F2,FCN4 

INTEGER  I, J, PASS, FLAG 
CHARACTER* 2  ANS 
EXTERNAL  FCN4 
C 

c -  DEFINE  CONSTANTS  - 

C 

SB  =  5.67051D-12 
K1  =  2 ,0D0*SB*EMIS 
K2  =  E*ABS 
C 

C -  SCALE  INPUTS  - 

C 

DELO  =  DEL0*100.0D0 
DELE  =  DELE* 1 00. ODO 
L  =  L*100.0D0 
K  =  K*0.01D0 
C 

C -  START  COUNTERS  - 

C 

PASS  =  0 
1  =  0 
J  =  1 
C 

c -  START  SEARCH  WITH  INTERVAL  [0.1, 0.2]  - 

C 

PRINT  200 
XI  =  O.lDO 
X2  =  0.2D0 

FI  =  FCN4(X1,TB,TE,HT, DELO, DELE, K,K1,K2,Q,L) 

5  1  =  1  +  1 
C 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  100  INTERVALS  - 

C 

IF(I.GT.J*100)  GOTO  10 

IF(PASS.EQ.O)  PRINT  100,1 

IF(PASS.EQ. 1)  PRINT  200 

IF(PASS.GE. 1)  PRINT  300 , XI . FI , X2 , F2 

F2  =  FCN4(X2,TB,TE,HT, DELO, DELE, K,K1,K2,Q,L) 

C 

C -  CHECK  FOR  FUNCTION  VALUES  OPPOSITE  IN  SIGN  - 

C 

IF(F1*F2.GE.0.0D0)  THEN 
XI  =  X2 
FI  =  F2 

IF(X2.LT.1.0D0)  X2  =  X2  +  O.lDO 

IF(X2.GE.1.0D0)  X2  =  X2  +  1 . ODO 
IF(PASS.GE. 1)  PASS  =  PASS  +  1 
GOTO  5 
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ENDIF 


C 

C - IF  INTERVAL  FOUND,  GO  FIND  ROOT - 

C 

PRINT  201 
GOTO  20 
C 

C - IF  NO  INTERVAL  FOUND,  CONTINl«:  SEARCH  ? - 

C 

10  PRINT  400,J*100 
READ  500, ANS 

IF(ANS.NE. • Y' .AND.ANS.NE. ’N* )  THEN 
PRINT  600 
GOTO  10 
ENDIF 

IF(ANS.EQ. 'N' )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 
C 

C -  IF  UNABLE  TO  FIND  INTERVAL,  LOOK  AT  FUNCTION  VALUES  ? 

C 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. 'Y' .AND.ANS.NE. ’N* )  THEN 
PRINT  600 
GOTO  15 
ENDIF 

IF{ANS.EQ. 'N* )  THEN 
FLAG  =  1 
RETURN 
ENDIF 
PASS  =  1 
1  =  0 
J  =  1 

XI  =  O.lDO 
X2  =  0.2D0 

FI  =  FCN4(Xl,TB,TE,HT,DEL0,DELE,K,Kl,K2,Q,L) 

GOTO  5 

C 

C -  GET  ROOT  IN  INTERVAL;  FIN  HEIGHT  AND  TEMP  AT  TIP  - 

C 

20  CALL  MDLIN4(FCN4,X1,X2,XR,TB,TE,HT,DEL0,DELE,K,K1,K2,Q,L) 
HT  =  XR 
C 

C -  CALCULATE  IDEAL  HEAT  DISSIPATED  - 

C 

QI  =  2.0DO*SB*EMIS*HT*L*TB**4.ODO 
C 
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c - CALCUIiATE  FIN  EFFICIENCY - 

C 

EFF  =  Q/QI 
C 

C -  SCALE  OUTPUTS  - 

C 

HT  =  HT*0.01D0 
DELO  =  DEL0*0.01D0 
DELE  =  DELE* 0.0 IDO 
L  =  L*0.01D0 
K  =  K*100.0D0 
C 

C -  FORMAT  - 

C 

100  FORMAT (•+•,• LOOKING  IN  INTERVAL:  *,13,’  ') 

200  FORMAT (/,•  COMMENCING  SEARCH  FOR  INTERVAL  TO  BRACKET  ROOT -  *,/) 

201  FORMAT (/,'  INTERVAL  FOUND,  WILL  NOW  LOOK  FOR  ROOT  -  ’,/) 

300  FORMAT('  XL  = ' , 2X, DlO . 4 , 2X, • F( XL)  = * , 2X,D10 . 4, 2X, 

&  'XR  =• ,2X,D10.4,2X, 'FIXR)  =‘,2X,D10.4) 

400  FORMAT(/,'  NO  ROOT  FOUND  IN  *,13, *  INTERVALS.  CONTINUE  SEARCH’  , 
&  •  (Y  OR  N)  ?•  ) 

500  FORMAT (A2) 

600  FORMAT(/,'  ****************  INCORRECT  ENTRY  I  *****************  ') 
700  FORMAT (/,’  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N)? 

&  •  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  IF 
&  •  F(XR)  IS  OF  THE  SAME  SIGN  AS  F(XL)  AND  ITS  ABSOLUTE  VALUE 
&  •  IS  INCREASING,  THERE  IS  LITTLE  PROBABILITY  OF  FINDING  AN  ’,/, 
&  •  INTERVAL.  THE  SAME  GOES  FOR  VALUES  OF  F(XR)  THAT  OSCILLATE’,/, 
&  ’  AROUND  A  CERTAIN  VALUE.  FOR  F(XR)  TO  BECOME  OPPOSITE  IN  ’,/, 
&  ’  SIGN  TO  F(XR),  IT  MUST  PASS  THROUGH  ZERO.  ') 

RETURN 
END 
C 

C5C******  TPSY;  FCN4(X,TB,TE,HT, DELO, DELE, K,K1,K2,Q,L)  *************** 

C 

C  PARAMETERS : 

C 

C  X  -  INDEPENDENT  VARIABLE  (BASE  DERIV) 

C  TB  -  BASE  TEMPERATURE 

C  TE  -  TIP  TEMPERATURE 

C  HT  -  FIN  HEIGHT 

C  DELO  -  FIN  THICKNESS  AT  BASE 

C  DELE  -  FIN  THICKNESS  AT  TIP 

C  K  -  FIN  CONDUCTIVITY 
C  K1  -  CONSTANT  =  2.0*SB*EHIS 
C  K2  -  CONSTANT  =  ABS*E 
C  Q  -  HEAT  INPUT  THRU  BASE 
C  L  -  FIN  LENGHT 
C 

c - 
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c 

REAL*8  FUNCTION  FCN4 (X, TB, TE, HT, DELO , DELE, K, K1 , K2 , Q, L ) 
REAL*8  X,X0,XFINAL,T0(2) ,TEND(2) ,F(2) ,H,TB,TE,HT,DELO,DELE, 


&  K,Kl,K2,Q,L,TSAVE,TOL 

EXTERNAL  DERIV4 
C 

C -  GUESS  HEIGHT  - 

C 

HT  =  X 
C 

C -  INITIALIZE  LEFT  END  POINT  BC*S:  FCN  AND  DERIVATIVE  KNOWN 

C 


T0(1)=  TB 

T0(2)=  -Q/(L*DEL0*K) 

C 

C -  INITIALIZE  RIGHT  END  POINT  BC*S  TO  0 

C 

TEND( 1)=0.0D0 
TEND(2)=0.0D0 
C 

C -  DEFINE  PARAMETERS  - 

C 

IF  (X.LT.l.ODO)  H  =  0.05D0 
IF  (X.GE.l.ODO)  H  =  O.lDO 
TOL  =  0.000 IDO 
XO  =  O.ODO 
XFINAL  =  HT 
TSAVE  =  1000. ODO 
C 

C -  USE  RUNGE  KUTTE  FELHBERG  TO  SHOOT  FROM  LEFT  END  BC’S  TO 

C  RIGHT  BC*S  AS  A  FUNCTION  OF  FIN  HEIGHT 

C 

10  IF(X0.LT.XFINAL+0.0001D0)  THEN 

CALL  RKFSY4 ( DER1V4 , XO , TO , TEND, F, H , HT, DELO , DELE , K, Kl , K2 , TOL ) 

C 

C - T  INCREASING  ME.ANS  DIVERGENCE - 

C 

IF  (TEND( 1) .GT. TSAVE)  GOTO  20 
T0(1)  =  TEND(l) 

TSAVE  =  TEND ( 1 ) 

TO (2)  =  TEND(2) 

GO  TO  10 
END  IF 

20  CONTINUE 

FCN4  =  TEND(2) 

TE  =  TEND(l) 

RETURN 

END 

C 

C5D  *****  TPSY:  RKFSY4 (DERIV4 , XO , TO , TEND, F, H, HT, DELO , DELE, K,K1 , K2 , TOL) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

1 


c 

c 

c 


PARAMETERS ; 

RKFSy4  -  SUBROUTINE  THAT  SOLVES  A  SYSTEM  OF  2  FIRST  ORDER 

DIFFERENTIAL  EQUATIONS  BY  THE  RUNGE-KUTTA-FEHLBERG 
METHOD.  THE  EQUATIONS  ARE  OF  THE  FORM: 

DT/DX  =  Y  =  Fl(X,T) 

DY/DX  =  F2(X,T,Y) 

DERIV4  -  A  SUBROUTINE  THAT  COMPUTES  VALUES  OF  THE  2  DERIVATIVES. 

XO  -  THE  INITIAL  VALUE  OF  INDEPENDENT  VARIABLE 
H  -  THE  INCREMENT  TO  T,  THE  STEP  SIZE 

TO  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 

TEND  -  AN  ARRAY  THAT  RETURNS  THE  FINAL  VALUES  OF  THE  FUNCTIONS 

F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 

H  -  STEP  SIZE 

HT  -  HEIGHT  OF  FIN 

DELO  -  THICKNESS  AT  BASE 

DELE  -  THICKNESS  AT  TIP 

K  -  FIN  CONDUCTIVITY 

Kl  -  CONSTANT  =  2 . 0*SB*EMIS 

K2  -  CONSTANT  =  E*ABS 

TOL  -  TOLERANCE 

TWRK  -  AN  ARRAY  USED  TO  HOLD  INTERMEDIATE  VALUES  DURING  THE 
COMPUTATION.  IT  MUST  BE  DIMENSIONED  OF  SIZE  6X2. 


SUBROUTINE  RKFSY4 ( DERIV4 , XO , TO , TEND , F , H , HT , DELO , DELE , K , Kl , K2 , TOL ) 
REAL*8  X0,T0(2) ,TEND(2) ,F(2) ,TWRK(6,2) , H,HT, DELO, DELE, K,Kl,K2, TOL, 
&  ERROR, SUM, STOREH,XEND 

INTEGER  I 

-  INITIALIZE  FOR  INTERVAL  [XO,  XO+H]  - 

XEND  =  XO+H 
STOREH  =  H 

-  CHECK  TO  SEE  IF  FINISHED  - 

IF(XO.GE.XEND)  THEN 
H  =  STOREH 
RETURN 
ELSE 


-  GET  FIRST  ESTIMATE  - 

CALL  DERIV4 (XO, TO, F,HT, DELO, DELE, K,K1,K2) 
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DO  10  I  =  1,2 

TWRK(1,I)  =  H*F(I) 

TEND(I)  =  T0(I)+TWRK(1,I)/4.0D0 
10  CONTINUE 
C 

C -  GET  SECOND  ESTIMATE  - 

C 

CALL  DERIV4(X0+H/4.0D0, TEND , F , HT , DELO , DELE , K , K 1 , K2 ) 

DO  20  I  =  1,2 

TWRK(2,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(TWRK(l,I)*3.OD0+TWRK(2,I)*9.0DO)/32.ODO 
20  CONTINUE 
C 

C -  GET  THIRD  ESTIMATE  - 

C 

CALL  DERIV4 (X0+3.0D0*H/8.0D0, TEND, F,HT, DELO, DELE, K,K1,K2) 

DO  30  I  =  1,2 

TWRK(3,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(TWRK(1,I)*1932.0D0-TWRK(2,I)*7200.0D0 
&  +  TWRK(3,I)*7296.0DO)/2197.0DO 

30  CONTINUE 
C 

c -  GET  FOURTH  ESTIMATE  - 

C 

CALL  DERIV4(X0+12 .0*H/ 13.0, TEND, F,HT, DELO, DELE, K,K1,K2) 

DO  40  I  =  1,2 

TWRK(4,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(439,0D0*TWRK(1,I>/216.0D0-8.0D0*TWRK(2,I) 
&  +  3680.0D0*TWRK(3,I)/513.0D0-845.0D0*TWRK(4,I)/4104.0D0) 

40  CONTINUE 
C 

c - get  fifth  estimate - . - 

c 

CALL  DERIV4 ( XO+H , TEND , F , HT , DELO , DELE , K , Kl , K2 ) 

DO  50  I  =  1,2 

TWRK(5,I)  =  H*F(I) 

TEND(I)  =  T0(I)-8.0D0*TWRK(1,I)/27.0+2.0D0*TWRK(2,I) 

&  -  3544.0D0*TWRK(3,I)/2565.0D0+1859.0D0*TWRK(4,I)/4104.0D0 

&  -  11.0D0*TWRK(5,I)/40.0D0 

50  CONTINUE 
C 

C -  GET  SIXTH  ESTIMATE  - 

C 

CALL  DERIV4 ( XO+H/2 . ODO , TEND , F , HT, DELO , DELE , K, Kl , K2 ) 

DO  60  I  =  1,2 

TWRK(6,I)  =  H*F(I) 


60  CONTINUE 
C 

C -  ESTIMATE  ERROR  BY  COMPUTING  DIFFERENCE  BETWEEN  FOURTH  AND 

C  FIFTH  ORDER  EQUATIONS 
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c 

SUM  =  O.ODO 
ERROR  =  O.ODO 
DO  70  I  =  1,2 

SUM  =  DABS(TWRK(1,I)/360.0DO-128.0DO*TWRK{3,I)/4275.0DO 

&  -  2197.0D0*TWRK(4,I)/75240.0D0+TWRK(5,I)/50.0D0 

&  +  2.0D0*TWRK(6,I)/55.0D0) 

ERROR  =  DMAXl (ERROR, SUM) 

70  CONTINUE 
C 

C -  IF  ERROR  LESS  THAN  TOLERANCE,  COMPUTE  X  AT  END  OF  - 

C  INTERVAL  FROM  A  WEIGHTED  AVERAGE  OF  SIX  ESTIMATES 

C 

IF ( ERROR. LT.TOL)  THEN 
DO  80  I  =  1,2 

TEND{I)  =  TO(I)+16.0DO*TWRK(1,I)/135.0DO 
&  +  6656.0D0*TWRK(3,1)/12825.0D0 

&  +  28561. 0D0*TWRK(4,I)/56430. 0-9. 0D0*TWRK(5,I)/50.0D0 

&  +  2.0DO*TWRK(6,I)/55.0DO 

TO (I)  =  TEND(I) 

80  CONTINUE 
XO  =  XO+H 
ENDIF 


C 

c -  IF  ERROR  GREATER  THAN  TOLERANCE,  REDUCE  STEP  AND  GO 

C  AGAIN 

C 


IF ( ERROR. GT.TOL)  THEN 
H  =  H/2.0D0 
ENDIF 
C 


C -  IF  ERROR  IS  SIGNIFICANLTY  LESS  THAN  TOLERANCE,  RELAX  STEP 

C 

IF ( ERROR. LT.H*TOL/l O.ODO)  THEN 
H  =  H*2.0D0 
ENDIF 
C 

C -  IF  OVERSHOOT,  REDUCE  STEP  - 

C 

IF(X0+H.GT.XEND)  THEN 
H  =  XEND-XO 


ENDIF 

C 

ENDIF 

C 

c 

GO  TO  1 
END 
C 

C5E******  TPSY:  DERIV4(X,T,F,HT,DEL0,DELE,K,Kl,K2)  ****************** 
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c 

C  PARAMETERS : 

C 

C  X  -  INDEPENDENT  VARIABLE 

C  T  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 

C  F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 

C  HT  -  HEIGHT  OF  FIN 

C  DELO  -  FIN  THICKNESS  AT  BASE 

C  DELE  -  FIN  THICKNESS  AT  TIP 

C  K  -  FIN  CONDUCTIVITY 

C  Kl  -  CONSTANT  =  2 . 0*SB*EMIS 

C  K2  -  CONSTANT  =  E*ABS 

C 

C - 

C 

SUBROUTINE  DERIV4 ( X , T , F , HT , DELO , DELE , K, Kl , K2 ) 

REAL*8  X,T(2) ,F(2) ,HT, DELO, DELE, K,K1,K2, DEL, DELP 
DEL(X)  =  DELO+(DELE-DELO)*X/HT 
DELP  =  (DELE-DELO) /HT 
F(l)  =  T(2) 

F(2)  =  -(DELP*T(2) )/DEL(X)+(Kl*T(l)**4.0D0-K2)/(K*DEL(X) ) 
RETURN 
END 
C 

C5F  *****  TPSY:  MDLIN4 ( FCN4 , Xl , X2 , XR, XTOL, FTOL,NLIM,TB, TE, HT, DELO , 
C  DELE,K,K1,K2,Q,L) 

C 

C  PARAMETERS : 

C 

C  FCN4  -  FUNCTION  THAT  COMPUTES  VALUES  FOR  F.  MUST  BE  DECLARED 
C  EXTERNAL  IN  CALLING  PROGRAM 

C  Xl,X2  -  INITIAL  VALUES  OF  X.  F(X)  MUST  CHANGE  SIGNS  AT  THESE 
C  POINTS 

C  XR  -  RETURNS  THE  ROOT  TO  THE  MAIN  PROGRAM 

C  XTOL  -  TOLERANCE  FOR  X 

C  FTOL  -  TOLERANCE  FOR  F 

C  NLIM  -  LIMIT  TO  NUMBER  OF  ITERATIONS 

C  TB  -  BASE  TEMPERATURE 

C  TE  -  TIP  TEMPERATURE 

C  HT  -  FIN  HEIGHT 

C  DELO  -  BASE  THICKNESS 

C  DELE  -  TIP  THICKNESS 

C  K  -  FIN  CONDUCTIVITY 

C  Kl  -  CONSTANT  =  2.0*SB*EMIS 

C  K2  -  CONSTANT  =  E*ABS 

C  Q  -  HEAT  INPUT  THRU  BASE 

C  L  -  FIN  LENGTH 

C 

C - 

C 
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SUBROUTINE  MDLIN4 ( FCN4 , Xl , X2 , XR,TB, TE, HT, DELO , DELE, K, Kl , K2 , Q,L) 
REAL*8  XERR, FSAVE, Fl , F2 , FR, Xl , X2 , XR, XTOL, FTOL, TB, TE, HT, DELO, 

&  DELE,K,K1,K2,Q,L,FCN4,X1SAV,X2SAV,F1SAV,F2SAV 

INTEGER  I, J, PASS 
CHARACTER* 2  ANS 

DATA  XTOL, FTOL/0. 000 IDO, O.OOOOIDO/ 

C 

C -  INITIALIZE  VALUES  — - 

C 

PRINT  200 
XISAV  =  Xl 
X2SAV  -  X2 

Fl  =  FCN4(X1,TB,TE,HT, DELO, DELE, K,K1,K2,Q,L) 

F2  =  FCN4(X2,TB,TE,HT, DELO, DELE, K,K1,K2,Q,L) 

FISAV  =  Fl 
F2SAV  =  F2 
FSAVE  =  Fl 
C 

C -  INITIALIZE  COUNTER  - 

C 

1  =  0 
J  =  1 
PASS  =  0 
C 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  50  ITERATIONS  - 

C 

5  1  =  1  +  1 

IF(I.GT.J*50)  GOTO  10 
XR  =  X2-F2*(X2-X1)/(F2-F1) 

FR  =  FCN4(XR,TB,TE,HT, DELO, DELE, K,K1,K2,Q,L) 

IF(PASS.EQ.O)  PRINT  100,1 
IF(PASS.EQ. 1)  PRINT  200 
IF(PASS.GE.l)  PRINT  300,I,XR,FR 
IF(PASS.GE. 1)  PASS  =  PASS  +  1 
C 

C -  CHECK  STOPPING  CRITERIA  - 

C 

XERR  =  DABS(X1-X2)/2.0D0 

IF(XERR.LE,XTOL.OR.DABS(5R) .LE.FTOL)  RETURN 
C 

C - FIND  NEW  POINT - 

C 

IF(FR*F1.GE,0.0D0)  THEN 
Xl  =  XR 
Fl  =  FR 

IF(FR*FSAVE.GT,0.0D0)  F2  =  F2/2.0D0 
FSAVE  =  FR 
ELSE 
X2  =  XR 
F2  =  FR 
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IF(FR*FSAVE.GT.O.ODO)  FI  =  F1/2.0D0 
FSAVE  =  FR 
ENDIF 
GOTO  5 
C 

C - FAIL  TO  FIND  ROOT,  CONTINUE  SEARCH  ? - 

C 

10  PRINT  400,J*50 
READ  500, ANS 

IF(ANS.NE. 'Y' .AND.ANS.NE. 'N* )  THEN 
PRINT  600 
GOTO  10 
ENDIF 

IF(ANS.EQ. 'N* )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 
C 

C -  FAIL  TO  FIND  ROOT,  LOOK  AT  FUNTION  VALUES  ? 

C 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. 'Y' .AND.ANS.NE. 'N' )  THEN 
PRINT  600 
GOTO  15 
ENDIF 

IF(ANS.EQ. 'N' )  STOP 
XI  =  XISAV 
X2  =  X2SAV 
Fl  =  FISAV 
F2  =  F2SAV 
1  =  0 
J  =  1 
PASS  =  1 
GOTO  5 
C 

C  -  FORMAT  - 

C 


100  FORMAT(’+',  ’ITERATION  :  ’,13,’  ') 

200  FORMAT (/,’  COMMENCING  SEARCH  FOR  ROOT  IN  INTERVAL....  ’,/) 

300  FORMAT(’  AT  ITERATION ', 13 , 3X, ’  X  =  •,D15.6,3X, '  F(X)  =  ’,015.6) 

400  FORMAT(/,’  FAILED  TO  FIND  ROOT  AFTER  ’,13,’  ITERATIONS.’  , 

&  ’  CONTINUE  SEARCH  (Y  OR  N)  ?’  ) 

500  FORMAT(A2) 

600  FORMAT (/’  *****************  INCORRECT  ANSWER  ********************  ) 
700  FORMAT (/,’  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N  )?  ’,/, 

&  ’  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  ’,/, 

&  ’  IF  F(XR)  IS  OSCILLATING  ABOUT  A  CERTAIN  VALUE  OR  ITS  ’,/, 

i  ’  ABSOLUTE  VALUE  IS  INCREASING  THEN  THERE  IS  LITTLE  ’,/, 

&  ’  PROBABILITY  THAT  A  ROOT  WILL  BE  FOUND.  ’) 
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END 


C 

C6A  *****  TRAN  HEADERS  *********************************************** 
C 

C  SUBROUTINE  TRAN (TB, TE , L, HT, DELO , DELE, K,ABS , EMIS , E, QI , Q, EFF, FLAG) 

C 

C - 

c 

C  FIN  SHAPE:  TRIANGULAR 
C 

C  PROBLEM  TYPE:  ANALYSIS 
C 

C  INPUT:  TB,L,HT,DELO,K,ABS,EMIS,E 
C 

C  OUTPUT :  TE , QI , Q , EFF 
C 

C - 

C 

C  PARAMETERS : 

C 

C  TB  -  TEMPERATURE  AT  THE  BASE  OF  THE  FIN 

C  TE  -  TEMPERATURE  AT  THE  TIP  OF  THE  FIN 

C  HT  -  HEIGHT  OF  THE  FIN 
CL-  LENGTH  OF  THE  FIN 
C  DELO  -  THICKNESS  AT  FIN  BASE 

C  DELE  -  THICKNESS  AT  FIN  TIP:  DELE  =  0.01*DELO 

C  K  -  CONDUCTIVITY  OF  THE  FIN 

C  ABS  -  ABSORPTIVITY  OF  THE  FIN 

C  EMIS  -  EMISSIVITY  OF  THE  FIN 

C  E  -  EXTERNAL  HEAT  INCIDENT  ON  THE  FIN 

C  QI  -  IDEAL  HEAT  DISSIPATED  BY  THE  FIN 

C  Q  -  REAL  HEAT  DISSIPATED  BY  THE  FIN 

C  EFF  -  EFFICIENCY  OF  THE  FIN 

C  FLAG  -  0  =  CONVERGENCE,  1  =  DIVERGENCE 

C 

C - 

c 

C  FUNCTIONS:  FCN5 ( X, TB,TE, L, HT, DELO , DELE, K, Kl ,K2 , Q) 

C 

C  FUNCTION  WHOSE  ROOT  IS  FOUND  (BASE  DERIV  &  TIP  TEMP) 

C 

c - 

c 

C  SUBROUTINES : 

C 

C  MDLIN4(FCN5,X1,X2,XR,XT0L,FT0L,NL1M,TB,HT,DEL,K,K1,K2)  : 

c 

C  FINDS  THE  ROOT  OF  AN  EQUATION  FCN4  BY  THE  METHOD  OF  MODIFIED 
C  LINEAR  INTERPOLATION 
C 


122 


C  RKFSY5(DERIV5,XBEGIN,H, TO, TEND, F,HT,DELO, DELE, K,K1,K2,T0L) : 

C 

C  SOLVES  SECOND  ORDER  NONLINEAR  DE  BY  RUNGE-KUTTE-FEHLBERG  METHOD 
C 

C  DERIV5(X,T,F,HT,DEL0,DELE,K,K1,K2) : 

C 

C  COMPUTES  DERIVATIVES  FOR  RKFSY5 
C 

C6B  *****  TRAN  MAINS  ************************************************* 

C 

SUBROUTINE  TRAN (TB, TE , L, HT, DELO, DELE, K,ABS, EMIS , E , QI , Q, EFF, FLAG) 
REAL*8  TB, TE,L,HT, DELO, DELE, K,ABS, EMIS, E,QI,Q, EFF, XI, X2,XR, 

&  SB,TBDER,K1,K2,F1,F2,FCN5 

INTEGER  I, J, PASS, FLAG 
CHARACTER* 2  ANS 
EXTERNAL  FCN5 
C 

c -  DEFINE  CONSTANTS  - 

C 

SB  =  5.67051D-12 
K1  =  2.0D0*SB*EMIS 
K2  =  E*ABS 
C 

C -  SCALE  INPUTS  - 

C 

HT  =  HT*100.0D0 
DELO  =  DEL0*100.0D0 
L  =  L*100.0D0 
K  =  K*0.01D0 
C 

C - FOR  TRIANGULAR  FIN:  TIP  THICKNESS  =  0.0 IDO  X  BASE  THICKNESS 

C 

DELE  =  DEL0*0.01D0 
C 

C -  START  COUNTERS  - 

C 

PASS  =  0 
1  =  0 
J  =  1 
c 

c -  START  SEARCH  WITH  INTERVAL  [-1,0]  - 

C 

PRINT  200 
XI  =  O.ODO 
X2  =  -l.ODO 

FI  =  FCN5(X1,TB,TE,HT, DELO, DELE, K,K1,K2) 

5  1  =  1  +  1 
C 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  100  INTERVALS  - 

C 
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IF(I.GT.J*100)  GOTO  10 
IF(PASS.EQ.O)  PRINT  100,1 
IF(PASS.EQ. 1)  PRINT  200 
IF(PASS.GE.l)  PRINT  300 , XI , Fl, X2,F2 


F2  =  FCN5(X2,TB,TE,HT,DEL0,DELE,K,K1,K2) 

C 

C -  CHECK  FOR  FUNCTION  VALUES  OPPOSITE  IN  SIGN 

C 


IF(F1*F2.GE.0.0D0)  then 
XI  =  X2 
Fl  =  F2 

X2  =  X2  -  l.ODO 
IF(PASS.GE.l)  PASS  »  PASS  +  1 
GOTO  5 
ENDIF 
C 

C - IF  INTERVAL  FOUND,  GO  FIND  ROOT - 

C 

PRINT  201 
GOTO  20 
C 

C -  IF  NO  INTERVAL  FOUND,  CONTINUE  SEARCH  ?  - 

C 

10  PRINT  400,J*100 
READ  500, ANS 

IF(ANS.NE. • y .AND.ANS-NE. 'N* )  THEN 
PRINT  600 
GOTO  10 
ENDIF 

IF(ANS.EQ. 'N* )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 
C 

C -  IF  UNABLE  TO  FIND  INTERVAL,  LOOK  AT  FUNCTION  VALUES  7 

C 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. ' Y' .AND.ANS.NE. ’N* )  THEN 
PRINT  600 
GOTO  15 
ENDIF 

IF(ANS.EQ. 'N' )  THEN 
FLAG  =  1 
RETURN 
ENDIF 
PASS  =  1 
1  =  0 
J  =  1 

XI  =  O.ODO 


124 


X2  =  -l.ODO 

FI  =  FCN5(X1,TB,TE,HT,DEL0,DELE,K,K1,K2) 

GOTO  5 
C 

C -  GET  ROOT  IN  INTERVAL:  DERIV  AT  BASE  AND  TEMP  AT  TIP 

C 

20  CALL  MDLIN5(FCN5,X1,X2,XR,TB,TE,HT,DEL0,DELE,K,K1,K2) 
TBDER  =  XR 


C 

C -  CALCULATE  IDEAL  HEAT  DISSIPATED  - 

C 

QI  =  2.0D0*SB*EMIS*HT*L*TB**4.ODO 
C 

C -  CALCULATE  REAL* 8  HEAT  DISSIPATED 

C 

Q  =  -K*DEL0*L*TBDER 
C 

C -  CALCULATE  FIN  EFFICIENCY  - 

C 


EFF  =  Q/QI 
C 
C 

C -  SCALE  OUTPUTS  - 

C 

DELO  =  DEL0*0.01D0 
DELE  =  DELE* 0.0 IDO 
L  =  L*0.01D0 
HT  =  HT*0.01D0 
K  =  K*100.0D0 
RETURN 
C 

c -  FORMAT  - 

C 

100  FORMAT( •+', 'LOOKING  IN  INTERVAL:  ',13,'  ') 

200  FORMAT (/,'  COMMENCING  SEARCH  FOR  INTERVAL  TO  BRACKET  ROCT _  ',/) 

201  FORMAT (/,'  INTERVAL  FOUND,  WILL  NOW  LOOK  FOR  ROOT  _  ',/) 

300  FORMATC  XL  = ' , 2X, DIO , 4 , 2X, ' F(XL)  = ' , 2X, DIO . 4 , 2X, 

&  'XR  =■ ,2X,D10.4,2X, •F(XR)  =',2X,D10.4) 

400  FORMAT (/,'  NO  ROOT  FOUND  IN  ',13,'  INTERVALS.  CONTINUE  SEARCH'  , 

5  '  (Y  OR  N)  ?'  ) 

500  FORMAT (A2) 

600  FORMAT(/,'  ****************  INCORRECT  ENTRY  I  *****************  ') 
700  FORMAT(/,'  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N)7  ',/, 

6  '  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  IF  ',/, 

&  '  F(XR)  IS  OF  THE  SAME  SIGN  AS  F(XL)  AND  ITS  ABSOLUTE  VALUE  ',/, 

&  '  IS  INCREASING,  THERE  IS  LITTLE  PROBABILITY  OF  FINDING  AN  ',/, 

&  '  INTERVAL.  THE  SAME  GOES  FOR  VALUES  OF  F(XR)  THAT  OSCILLATE',/, 

&  '  AROUND  A  CERTAIN  VALUE.  FOR  F(XR)  TO  BECOME  OPPOSITE  IN  ',/, 

&  '  SIGN  TO  F(XR),  IT  MUST  PASS  THROUGH  ZERO.  ') 
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END 


C 

C6C***»**  TRAN:  FCN5(X,TB,TE,HT,DEL0,DELE,K,K1,K2)  ****************** 
C 

C  PARAMETERS : 

C 

C  X  -  INDEPENDENT  VARIABLE  (BASE  DERIV) 

C  TB  -  BASE  TEMPERATURE 

C  TE  -  TIP  TEMPERATURE 

C  HT  -  FIN  HEIGHT 

C  DELO  -  FIN  THICKNESS  AT  BASE 

C  DELE  -  FIN  THICKNESS  AT  TIP 

C  K  -  FIN  CONDUCTIVITY 
C  K1  -  CONSTANT  =  2.0*SB*EMIS 
C  K2  -  CONSTANT  =  ABS*E 
C 

c - 

c 

REAL*8  FUNCTION  FCN5 ( X, TB , TE , HT, DELO , DELE, K, Kl , K2 ) 

REAL* 8  X,X0,XFINAL,T0(2) , TEND( 2 ) , F ( 2 ) , H , TB, TE , HT, DELO , DELE, 


&  K,K1,K2,TSAVE,T0L 

EXTERNAL  DERIV5 
C 

C -  INITIALIZE  LEFT  END  POINT  BC’S:  KNOW  FCT  VALUE,  GUESS 

C  DERIVATIVE 


C 

T0( 1)=  TB 
T0(2)=  X 
C 

C - INITIALIZE  RIGHT  END  POINT  BC’S  TO  0  - 

C 

TEND( 1)=0.0D0 
TEND(2)=0.0D0 
C 

C -  DEFINE  PARAMETERS  - 

C 

H=  O.lDO 
TOL  =  O.OOOIDO 
XO  =  O.ODO 
XFINAL  =  HT 
TSAVE  =  1000. ODO 
C 

C -  USE  RUNGE  KUTTE  FELHBERG  TO  SHOOT  FROM  LEFT  END  BC’S  TO 

C  RIGHT  BC’S  AS  A  FUNCTION  OF  BASE  DERIV 

C 

10  IF(X0.LT.XFINAL+0.0001D0)  THEN 

CALL  RKFSYS (DERIV5 , XO , TO , TEND, F, H, HT, DELO , DELE, K, Kl , K2 , TOL ) 

C 

C -  INCREASING  T  MEANS  DIVERGENCE  - 

C 
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20 


C 

C6D 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c — 
c 


c 

c 

c 


IF  (TEND( 1) .GE.TSAVE)  GOTO  20 
T0(1)  =  TEND(l) 

TSAVE  =  TEND(l) 

TO (2)  =  TEND(2) 

GO  TO  10 
ENDIF 
CONTINUE 
FCN5  =  TEND (2) 

TE  =  TEND(l) 

RETURN 

END 


*****  TRAN;  RKFSY5(DERIV5,X0, TO, TEND, F,H,HT,DEL0, DELE, K,Kl,K2,TOL) 
PARAMETERS : 

RKFSY5  -  SUBROUTINE  THAT  SOLVES  A  SYSTEM  OF  2  FIRST  ORDER 

DIFFERENTIAL  EQUATIONS  BY  THE  RUNGE-KUTTA-FEHLBERG 
METHOD.  THE  EQUATIONS  ARE  OF  THE  FORM: 

DT/DX  =  Y  =  F1(X,T) 

DY/DX  =  F2(X,T,Y) 

DERIV5  -  A  SUBROUTINE  THAT  COMPUTES  VALUES  OF  THE  2  DERIVATIVES. 

XO  -  THE  INITIAL  VALUE  OF  INDEPENDENT  VARIABLE 

TO  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 

TEND  -  AN  ARRAY  THAT  RETURNS  THE  FINAL  VALUES  OF  THE  FUNCTIONS 

F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 

H  -  STEP  SIZE 

HT  -  HEIGHT  OF  FIN 

DELO  -  THICKNESS  AT  BASE 

DELE  -  THICKNESS  AT  TIP 

K  -  FIN  CONDUCTIVITY 

K1  -  CONSTANT  =  2.0*SB*EMIS 

K2  -  CONSTANT  =  E*ABS 

TOL  -  TOLERANCE 

TWRK  -  AN  ARRAY  USED  TO  HOLD  INTERMEDIATE  VALUES  DURING  THE 
COMPUTATION.  IT  MUST  BE  DIMENSIONED  OF  SIZE  6X2. 


SUBROUTINE  RKFSY5 ( DERIV5 , XO , TO , TEND, F, H , HT , DELO , DELE , K, K1 , K2 , TOL ) 
REAL* 8  X0,T0(2) , TEND( 2 ) , F( 2 ) ,TWRK( 6 , 2 ) , H, HT,DEL0 , DELE, K,Kl , K2 , 

&  TOL, ERROR, SUM, STOREH,XEND 

INTEGER  I 

-  INITIALIZE  FOR  INTERVAL  [XO,  XO+H]  - 


XEND  =  XO+H 
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STOREH  =  H 


C 

C - CHECK  TO  SEE  IF  WE  ARE  FINISHED 

C 

1  IF(XO.GE,XEND)  THEN 


H  =  STOREH 
RETURN 
ELSE 
C 

c - get  first  estimate  of  the  delta  X‘S - 

c 

CALL  DERIV5 ( XO , TO , F, HT, DELO ,DELE, K, K1 , K2 ) 

DO  10  I  =  1,2 

TWRK(1,I)  =  H*F(I) 

TEND(I)  =  T0(I)+TWRK(1,I)/4.0D0 
10  CONTINUE 
C 

c -  get  second  estimate  - 

c 

CALL  DERIV5 ( XO+H/4 . ODO , TEND, F, HT, DELO , DELE , K, K1 , K2 ) 

DO  20  I  =  1,2 

TWRK(2,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(TWRK(1,I)*3.0D0+TWRK(2,I)*9.0D0)/32.0D0 
20  CONTINUE 
C 

c -  get  third  estimate  - 

c 

CALL  DERIV5 ( XO+3 . 0D0*H/ 8 . ODO , TEND , F, HT, DELO , DELE, K , K1 , K2 ) 

DO  30  I  =  1,2 

TWRK(3,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(TWRK(1,I)*1932.0D0-TWRK(2,I)*7200.0D0 
&  +  TWRK(3,I)*7296.ODO)/2197.OD0 

30  CONTINUE 
C 

c -  get  fourth  estimate  - 

c 

CALL  DERIV5 (X0+12,0*H/13.0, TEND , F, HT, DELO , DELE , K, Kl , K2 ) 

DO  40  I  =  1,2 

TWRK(4,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(439.0D0*TWRK(1,I)/216,0D0-8.0D0*TWRK(2,I) 


&  +  3680. 0D0*TWRK(3,I)/513.0D0-845.0D0*TWRK(4,I)/4104. ODO) 

40  CONTINUE 
C 

c -  get  fifth  estimate  - 


c 

CALL  DERIV5 ( XO+H , TEND , F , HT , DELO , DELE , K, Kl , K2 ) 

DO  50  I  =  1,2 

TWRK(5,I)  =  H*F(I) 

TEND(I)  =  T0(I)-8.0D0*TWRK(1,I)/27.0+2.0D0*TWRK(2,I) 

&  -  3544. ODO*TWRK(3,I)/2565.0DO+1859.0DO*TWRK(4,I)/4104. ODO 


128 


& 


-  11.0D0*TWRK(5,I)/40.0D0 
50  CONTINUE 
C 

C -  GET  SIXTH  ESTIMATE  - 


C 

CALL  DERIV5 ( XO+H/2 . ODO , TEND , F , HT , DELO , DELE , K , K1 , K2 ) 

DO  60  I  =  1,2 

TWRK(6,I)  =  H*F(I) 

6  0  CONTINUE 

C 

C -  ESTIMATE  THE  ERROR  BY  COMPUTING  THE  DIFFERENCE  BETWEEN 

C  THE  FOURTH  AND  FIFTH  ORDER  EQUATIONS 

C 

SUM  =  O.ODO 
ERROR  =  O.ODO 
DO  70  I  =  1,2 

SUM  =  DABS(TWRK(1,I)/360.0D0-128.0D0*TWRK(3,I)/4275.0D0 

&  -  2197. 0D0*TWRK(4,I)/75240.0D0+TWRK(5,I)/50. ODO 

&  +  2.0D0*TWRK^6,I)/55.0D0) 

ERROR  =  DMAXl (ERROR, SUM) 

70  CONTINUE 

C 

C -  IF  ERROR  LESS  THAN  TOLERANCE,  COMPUTE  X  AT  END  OF  - 

C  THE  INTERVAL  FROM  A  WEIGHTED  AVERAGE  OF  SIX  ESTIMATES, 

C  THEN  RETURN . 

C 

IF ( ERROR. LT.TOL)  THEN 
DO  80  I  =  1,2 

TEND(I)  =  TO(I)+16.0DO*TWRK(1,I)/135.0DO 
&  +  6656. 0D0*TWRK(3,I)/12825. ODO 

&  +  28561. 0D0*TWRK(4,I)/56430. 0-9. 0D0*TWRK(5,I)/50. ODO 

&  +  2.0D0*TWRK(6,I)/55.0D0 

T0(I)  =  TEND(I) 

80  CONTINUE 
XO  =  XO+H 
END  IF 
C 

C -  IF  ERROR  GREATER  THAN  TOLERANCE,  THEN  REDUCE  STEP  - 

C 

IF < ERROR. GT.TOL)  THEN 
H  =  H/2.0D0 
END  IF 
C 


C -  IF  ERROR  IS  SIGNIFICANTLY  LESS  THAN  TOLERANCE,  RELAX  STEP 

C 

IF(ERROR.LT.H*TOL/10.0DO)  then 
H  =  H*2.0D0 
END  IF 
C 

C -  IF  OVERSHOOT,  REDUCE  STEP  - 
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c 

IF(XO+H.GT.XEND)  THEN 
H  =  XEND-XO 
ENDIF 
C 

ENDIF 

C 

GO  TO  1 
END 
C 

c 

C6E******  TRAN:  DERIV5 (X,T, F, HT,DEL0,DELE,K,K1 ,K2 )  ****************** 
C 

C  PARAMETERS : 

C 

C  X  -  INDEPENDENT  VARIABLE 

C  T  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 

C  F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 

C  HT  -  HEIGHT  OF  FIN 

C  DELO  -  FIN  THICKNESS  AT  BASE 

C  DELE  -  FIN  THICKNESS  AT  TIP 

C  K  -  FIN  CONDUCTIVITY 

C  K1  -  CONSTANT  =  2 . 0D0*SB*EMIS 

C  K2  -  CONSTANT  =  E*DABS 

C 

C - 

c 

S UBROUTI NE  DERI V5 ( X , T , F , HT , DEL 0 , DELE , K , K 1 , K2 ) 

REAL*8  X,T(2) ,F(2) ,HT, DELO, DELE, K,K1,K2, DEL, DELP 
DEL(X)  =  DEL0+( DELE -DELO )*X/HT 
DELP  =  (DELE -DELO )/HT 
F(l)  =  T(2) 

F(2)  =  -(DELP*T(2) )/DEL(X)+(K1*T(1)**4.0D0-K2)/(K*DEL(X) ) 

RETURN 

END 

C 

C6F  *****  TRAN:  MDLIN5(FCN5,Xl,X2,XR,XTOL,FTOL,NLIM,TB,TE,HT,DELO, 

C  DELE,K,K1,K2) 

C 

C  PARAMETERS : 

C 

C  FCN5  -  FUNCTION  THAT  COMPUTES  VALUES  FOR  F.  MUST  BE  DECLARED 
C  EXTERNAL  IN  CALLING  PROGRAM 

C  X1,X2  -  INITIAL  VALUES  OF  X.  F(X)  MUST  CHANGE  SIGNS  AT  THESE 
C  POINTS 

C  XR  -  RETURNS  THE  ROOT  TO  THE  MAIN  PROGRAM 
C  XTOL  -  TOLERANCE  FOR  X 
C  FTOL  -  TOLERANCE  FOR  F 
C  NLIM  -  LIMIT  TO  NUMBER  OF  ITERATIONS 
C  TB  -  BASE  TEMPERATURE 
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C  TE  -  TIP  TEMPERATURE 
C  HT  -  FIN  HEIGHT 
C  DELO  -  BASE  THICKNESS 
C  DELE  -  TIP  THICKNESS 
C  K  -  FIN  CONDUCTIVITY 
C  K1  -  CONSTANT  =  2.0*SB*EMIS 
C  K2  -  CONSTANT  =  E*ABS 
C 

C - 

c 

SUBROUTINE  MDLIN5 { FCN5 , Xl , X2 , XR, TB , TE , HT , DELO , DELE , K , Kl , K2 ) 
REAL* 8  XERR, FSAVE , F 1 , F2 , FR, Xl , X2 , XR, XTOL , FTOL , TB , TE , HT , DELO , 
&  DELE , K, Kl , K2 , FCN5 , XlSAV, X2SAV, FISAV, F2SAV 

INTEGER  I, J, PASS 
CHARACTER* 2  ANS 

DATA  XTOL,FTOL/0.0001D0, O.OOOOIDO/ 

C 

C -  INITIALIZE  VALUES  - 

C 

PRINT  200 
XlSAV  =  Xl 
X2SAV  =  X2 

FI  =  FCN5(X1,TB,TE,HT, DELO, DELE, K,K1,K2) 

F2  =  FCN5(X2,TB,TE,HT, DELO, DELE, K,K1,K2) 

FISAV  =  FI 
F2SAV  =  F2 
FSAVE  =  FI 
C 

C -  INITIALIZE  COUNTER  - 

C 

1  =  0 
J  =  1 
PASS  =  0 
C 

c -  LIMIT  SEARCH  TO  INCREMENTS  OF  50  ITERATIONS  - 

C 

5  1  =  1  +  1 

IF(I.GT.J*50)  GOTO  10 
XR  =  X2-F2*(X2-X1)/(F2-F1) 

FR  =  FCN5(XR,TB,TE,HT, DELO, DELE, K,K1,K2) 

IF(PASS.EQ.O)  PRINT  100,1 
IF(PASS.EQ. 1)  PRINT  200 
IF(PASS.GE. 1)  PRINT  300,I,XR,FR 
IF(PASS.GE. 1)  PASS  =  PASS  +  1 
C 

C -  CHECK  STOPPING  CRITERIA  - 

C 

XERR  =  DABS(X1-X2)/2.0D0 

IF(XERR.LE.XTOL.OR.DABS(FR) .LE.FTOL)  RETURN 
C 
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c - FIND  NEW  POINT - 

C 

IF(FR*F1.GE.0.0D0)  THEN 
XI  =  XR 
Fl  =  FR 

IF  (FR*FSAVE.GT.0.0D0)  F2  =  F2/2.0DO 
FSAVE  =  FR 
ELSE 
X2  =  XR 
F2  =  FR 

IF  (FR*FSAVE.GT.0.0D0)  Fl  =  F1/2.0D0 
FSAVE  =  FR 
ENDIF 
GOTO  5 
C 

C - FAIL  TO  FIND  ROOT,  CONTINUE  SEARCH  ? - 

C 

10  PRINT  400,J*50 
READ  500, ANS 

IF(ANS.NE. ’Y* .AND.ANS.NE. ‘N* )  THEN 
PRINT  600 
GOTO  10 
ENDIF 

IF(ANS.EQ. 'N' )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 
C 

c -  fail  to  find  root,  look  at  funtion  values  7  - 

c 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. 'Y' .AND.ANS.NE. ‘N’ )  THEN 
PRINT  600 
GOTO  15 
ENDIF 

IF(ANS.EQ. 'N' )  STOP 
XI  =  XlSAV 
X2  =  X2SAV 
Fl  =  FlSAV 
F2  =  F2SAV 
1  =  0 
J  =  1 
PASS  =  1 
GOTO  5 
C 

C  - 

c 

100  FORMAT('+*,  ’ITERATION  :  ’,13,’  •) 

200  FORMAT (/,'  COMMENCING  SEARCH  FOR  ROOT  IN  INTERVAL _  ’,/) 
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300  FORMATC  AT  ITERATION 13, 3X, •  X  =  •,D15.6,3X, *  F(X)  =  ',015.6) 


400  FORMAT(/,'  FAILED  TO  FIND  ROOT  AFTER  ',13,'  ITERATIONS,'  , 

&  '  CONTINUE  SEARCH  (Y  OR  N)  ?'  ) 

500  FORMAT (A2) 

600  FORMAT(/'  *****************  INCORRECT  ANSWER  *******************•) 
700  FORMAT (/,'  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N  )7  ',/, 

&  '  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  ',/, 

&  '  IF  F(XR)  IS  OSCILLATING  ABOUT  A  CERTAIN  VALUE  OR  ITS  ',/, 

&  '  ABSOLUTE  VALUE  IS  INCREASING  THEN  THERE  IS  LITTLE  ',/, 

&  '  PROBABILITY  THAT  A  ROOT  WILL  BE  FOUND.  ') 


END 

C 

C7A  *****  TRSY  HEADER6  *********************************************** 
C 

C  SUBROUTINE  TRSY ( TB , TE , HT , L , DELO , DELE, K, DABS , EMIS , E , QI , 0, EFF, FLAG ) 

C 

C - 

c 

C  FIN  SHAPE:  TRIANGULAR 
C 

C  PROBLEM  TYPE;  SYNTHESIS 
C 

C  INPUT;  TB,L, DELO, K,ABS, EMIS, E,Q 
C 

C  OUTPUT;  TE,HT,QI,EFF 
C 

C - 

c 

C  PARAMETERS : 

C 

C  TB  -  TEMPERATURE  AT  THE  BASE  OF  THE  FIN 

C  TE  -  TEMPERATURE  AT  THE  TIP  OF  THE  FIN 

C  HT  -  HEIGHT  OF  THE  FIN 
CL-  LENGTH  OF  THE  FIN 
C  DELO  -  THICKNESS  AT  FIN  BASE 

C  DELE  -  THICKNESS  AT  FIN  TIP;  DELE  =  0.01D0*DEL0 

C  K  -  CONDUCTIVITY  OF  THE  FIN 

C  DABS  -  ABSORPTIVITY  OF  THE  FIN 

C  EMIS  -  EMISSIVITY  OF  THE  FIN 

C  E  -  EXTERNAL  HEAT  INCIDENT  ON  THE  FIN 

C  QI  -  IDEAL  HEAT  DISSIPATED  BY  THE  FIN 

C  Q  -  REAL  HEAT  DISSIPATED  BY  THE  FIN 

C  EFF  -  EFFICIENCY  OF  THE  FIN 

C  FLAG  -  0  =  CONVERGENCE,  1  =  DIVERGENCE 

C 

C - 

c 

C  FUNCTIONS;  FCN6 (X,TB,TE, L,HT, DELO, DELE, K,Kl ,K2 , Q) 

C 
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C  FUNCTION  WHOSE  ROOT  IS  FOUND  (BASE  DERIV  &  TIP  TEMP) 

C 

C - 

c 

C  SUBROUTINES ; 

C 

C  MDLIN6 ( FCN6 , Xl , X2 , XR, XTOL , FTOL , NLIM, TB , HT , DEL , K , Kl , K2 ) : 

C 

C  FINDS  THE  ROOT  OF  AN  EQUATION  FCN4  BY  THE  METHOD  OF  MODIFIED 
C  LINEAR  INTERPOLATION 
C 

C  RKFSY6 ( DERIV6 , XBEGIN , H , TO , TEND, F, HT, DELO , DELE , K, Kl , K2 , TOL ) ; 

C 

C  SOLVES  SECOND  ORDER  NONLINEAR  DE  BY  RUNGE-KUTTE-FEHLBERG  METHOD 
C 

C  DERIV6 (X,T,F,HT,DEL0,DELE,K,K1,K2) ; 

C 

C  COMPUTES  DERIVATIVES  FOR  RKFSY6 
C 

C7B  *****  TRSY  MAIN6  ************************************************* 
C 

SUBROUTINE  TRSY ( TB , TE , L, HT, DELO , DELE, K,ABS , EMIS , E , QI , Q, EFF, FLAG ) 
REAL*8  TB, TE,L,HT, DELO, DELE, K,ABS, EMIS, E,QI,Q, EFF, Xl,X2,XR, 

&  SB,K1,K2,F1,F2,FCN6 

INTEGER  I, J, PASS, FLAG 
CHARACTER*2  ANS 
EXTERNAL  FCN6 
C 

C -  DEFINE  CONSTANTS  - 

C 

SB  =  5.67051D-12 
Kl  =  2.0D0*SB*EMIS 
K2  =  E*ABS 
C 

C -  SCALE  INPUTS  - 

C 

DELO  =  DEL0*100.0D0 
L  =  L*100.0D0 
K  =  K*0.01D0 
C 

C -  FOR  TRIANGULAR  FIN;  TIP  THICKNESS  =  0.0 IDO  BASE  THICKNESS  — 

C 

DELE  =  DELO *0.0 IDO 
C 

C -  START  COUNTERS  - 

C 

PASS  =  0 
1  =  0 
J  =  1 
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c 

C -  START  SEARCH  WITH  INTERVAL  [0.1, 0.2]  - 

C 

PRINT  200 
XI  =  l.ODO 
X2  =  2.0D0 

Fl  =  FCN6 (X1,TB,TE,HT,DEL0,DELE,K,K1,K2,Q,L) 

5  1  =  1  +  1 
C 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  100  INTERVALS 


C 

IF(I.GT.J*100)  GOTO  10 

IF(PASS.EQ.O)  PRINT  100,1 

IF(PASS.EQ. 1)  PRINT  200 

IF(PASS.GE. 1)  PRINT  300 , XI , Fl , X2 , F2 

F2  =  FCN6(X2,TB,TE,HT,DEL0,DELE,K,K1,K2,Q,L) 

C 

C -  CHECK  FOR  FUNCTION  VALUES  OPPOSITE  IN  SIGN  - 

C 

IF(F1*F2.GE.0.0D0)  THEN 
XI  =  X2 
Fl  =  F2 

X2  =  X2  +  l.ODO 

IF(PASS.GE.l)  PASS  =  PASS  +  1 
GOTO  5 
ENDIF 
C 

c - IF  INTERVAL  FOUND,  GO  FIND  ROOT - 

C 

PRINT  201 
GOTO  20 
C 

C -  IF  NO  INTERVAL  FOUND,  CONTINUE  SEARCH  ?  - 

C 

10  PRINT  400,J*100 
READ  500, ANS 

IF(ANS,NE. ’Y* .AND.ANS.NE. ’N’ )  THEN 
PRINT  600 
GOTO  10 
ENDIF 

IF(ANS.EQ. 'N' )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 
C 

r -  IF  UNABLE  TO  FIND  INTERVAL,  LOOK  AT  FUNCTION  VALUES  ? 

C 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. • Y’ .AND.ANS.NE. ‘N’ )  THEN 
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PRINT  600 
GOTO  15 
ENDIF 

IF(ANS.EQ. 'N* )  THEN 
FLAG  =  1 
RETURN 
ENDIF 
PASS  =  1 
1  =  0 
J  =  1 

XI  =  l.ODO 
X2  =  2.0D0 

FI  =  FCN6(X1,TB,TE,HT,DEL0,DELE,K,K1,K2,Q,L) 

GOTO  5 
C 

C -  GET  ROOT  IN  INTERVAL:  FIN  HEIGHT  AND  TEMP  AT  TIP  - 

C 

20  CALL  MDLIN6(FCN6,X1,X2,XR,TB,TE,HT,DEL0,DELE,K,K1,K2,Q,L) 

HT  =  XR 
C 

C -  CALCULATE  IDEAL  HEAT  DISSIPATED  - 

C 

QI  =  2.0D0*SB*EMIS*HT*L*TB»»4.0D0 

c 

C -  CALCULATE  FIN  EFFICIENCY  - 

C 

EFF  =  0/QI 
C 

C -  SCALE  OUTPUTS  - 

C 

HT  =  HT*0.01D0 
DELO  =  DEL0*0.01D0 
L  =  L*0.01D0 
K  =  K*100.0D0 
RETURN 
C 

c -  FORMAT  - 

C 

100  FORMAT( '+•, 'LOOKING  IN  INTERVAL:  ',13,'  ') 

200  FORMAT (/,'  COMMENCING  SEARCH  FOR  INTERVAL  TO  BRACKET  ROOT _  ',/) 

201  FORMAT (/,'  INTERVAL  FOUND,  WILL  NOW  LOOK  FOR  ROOT  _  ',/) 

300  FORMATC  XL  = ' , 2X, DIO . 4 , 2X, •F(XL)  =' ,2X,D10.4,2X, 

&  'XR  =• ,2X,D10.4,2X, •F(XR)  =',2X,D10.4) 

400  FORMAT(/,'  NO  ROOT  FOUND  IN  ',13,'  INTERVALS.  CONTINUE  SEARCH'  , 
i  '  (Y  OR  N)  ?•  ) 

500  F0RMAT(A2) 

600  FORMAT(/,'  ****************  INCORRECT  ENTRY  I  *****************  •) 
700  FORMAT (/,'  WOULD  YOU  LIKE  TO  SEE  FUNCTION  VALUES  (Y  OR  N)? 

&  •  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  IF  ',/, 
&  '  F(XR)  IS  OF  THE  SAME  SIGN  AS  F(XL)  AND  ITS  ABSOLUTE  VALUE  ',/, 
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&  •  IS  INCREASING,  THERE  IS  LITTLE  PROBABILITY  OF  FINDING  AN 
&  •  INTERVAL.  THE  SAME  GOES  FOR  VALUES  OF  F<XR)  THAT  OSCILLATE',/, 
&  •  AROUND  A  CERTAIN  VALUE.  FOR  F(XR)  TO  BECOME  OPPOSITE  IN  ’,/, 
&  •  SIGN  TO  F(XR),  IT  MUST  PASS  THROUGH  ZERO,  ') 

END 
C 

C7C******  TRSY;  FCN6 (X, TB,TE, HT,DEL0 ,DELE,K, Kl,K2 , Q,L)  ************** 

C 

C  PARAMETERS : 

C 

C  X  -  INDEPENDENT  VARIABLE  (BASE  DERIV) 


c 

TB  - 

BASE  TEMPERATURE 

c 

TE  - 

TIP  TEMPERATURE 

c 

HT  - 

FIN  HEIGHT 

c 

DELO 

-  FIN  THICKNESS  AT  BASE 

c 

DELE 

-  FIN  THICKNESS  AT  TIP 

c 

K  - 

FIN  CONDUCTIVITY 

c 

Kl  - 

CONSTANT  =  2.0*SB*EMIS 

c 

K2  - 

CONSTANT  =  ABS*E 

c 

Q  - 

HEAT  INPUT  THRU  BASE 

c 

L  - 

FIN  LENGTH 

C 

c - 

c 

REAL*8  FUNCTION  FCN6 ( X, TB , TE, HT, DELO , DELE , K, Kl , K2 , Q, L) 
REAL*8  X,X0,XFINAL,T0(2) ,TEND(2) ,F(2) ,H,TB,TE,HT,DEL0,DELE, 


&  K,Kl,K2,Q,L,TSAVE,TOL 

EXTERNAL  DERIV6 
C 

C -  GUESS  HEIGHT  - 

C 

HT  =  X 
C 

c -  INITIALIZE  LEFT  END  POINT  BC'S:  FCN  AND  DERIVATIVE  KNOWN 

C 


T0( 1)=  TB 

T0(2)=  -Q/(L*DEL0*K) 

C 

c -  INITIALIZE  RIGHT  END  POINT  BC'S  TO  0 

C 

TEND( 1)=0.0D0 
TEND(2)=0.0D0 
C 


c -  DEFINE  PARAMETERS  - 

C 

C  IF  (X.LT.l.ODO)  H  =  O.OIDO 

C  IF  (X.GE.l.ODO)  H  =  0 , IDO 

H  =  O.OIDO 
TOL  =  O.OOOIDO 
XO  =  O.ODO 
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XFINAL  =  X 
TSAVE  =  1000. ODO 
C 

C -  USE  RUNGE  KUTTE  FELHBERG  TO  SHOOT  FROM  LEFT  END  BC • S  TO 

C  RIGHT  BC'S  AS  A  FUNCTION  OF  FIN  HEIGHT 

C 

10  IF(X0.LT.XFINAL+0.0001D0)  THEN 

CALL  RKFSVe ( DERIV6 , XO , TO , TEND, F, H, HT, DELO , DELE , F, Kl , K2 , TOL) 

C 

C -  T  INCREASING  MEANS  DIVERGENCE  - 

C 

IF  (TEND( 1) .GT. TSAVE)  GOTO  20 
T0(1)  =  TEND(l) 

TSAVE  =  TEND(l) 

TO (2)  =  TEND (2) 

GO  TO  10 
ENDIF 

20  CONTINUE 

FCN6  =  TEND(2) 

TE  =  TEND { 1 ) 

RETURN 

END 

C 

C7D  *****  TRSYs  RKFSY6 (DERIV6,X0, TO, TEND, F,H,HT, DELO, DELE, K,K1,K2, TOL) 
C 

C  PARAMETERS  s 
C 

C  RKFSY6  -  SUBROUTINE  THAT  SOLVES  A  SYSTEM  OF  2  FIRST  ORDER 
C  DIFFERENTIAL  EQUATIONS  BY  THE  RUNGE -KUTTA-FEHLBERG 

C  METHOD.  THE  EQUATIONS  ARE  OF  THE  FORM: 

C 

C  DT/DX  =  Y  =  Fl(X,T) 

C  DY/DX  =  F2(X,T,Y) 

C 

C  DERIV6  -  A  SUBROUTINE  THAT  COMPUTES  VALUES  OF  THE  2  DERIVATIVES. 

C 

C  XO  -  THE  INITIAL  VALUE  OF  INDEPENDENT  VARIABLE 

C  TO  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 
C  TEND  -  AN  ARRAY  THAT  RETURNS  THE  FINAL  VALUES  OF  THE  FUNCTIONS 
C  F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 
C  H  -  STEP  SIZE 
C  HT  -  HEIGHT  OF  FIN 
C  DELO  -  THICKNESS  AT  BASE 
C  DELE  -  THICKNESS  AT  TIP 
C  K  -  FIN  CONDUCTIVITY 
C  K1  -  CONSTANT  =  2.0*SB*EMIS 
C  K2  -  CONSTANT  =  E*ABS 
C  TOL  -  TOLERANCE 

C  TWRK  -  AN  ARRAY  USED  TO  HOLD  INTERMEDIATE  VALUES  DURING  THE 
C  COMPUTATION.  IT  MUST  BE  DIMENSIONED  OF  SIZE  6X2. 
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c 

c - 

c 

SUBROUTINE  RKFSY 6 ( DERIV6 , XO , TO , TEND , F , H , HT , DELO , DELE , K , K 1 , K2 ,  TOL ) 
REAL*8  X0,T0(2) ,TEND(2) ,F(2) ,TWRK(6,2) ,H,HT,DELO,DELE,K,K1,K2,TOL, 


&  ERROR, SUM, STOREH,XEND 

INTEGER  I 
C 

C - INITIALIZE  FOR  INTERVAL  [XO,  XO+H] 

C 


XEND  =  XO+H 
STOREH  =  H 


C 

C - CHECK  TO  SEE  IF  FINISHED 

C 

1  IF(XO.GE.XEND)  THEN 


H  =  STOREH 
RETURN 
ELSE 
C 

c -  get  first  ESTIMATE  - 

C 

CALL  DERIVC (XO, TO, F,HT, DELO, DELE, K,K1,K2) 

DO  10  I  =  1,2 

TWRK(1,I)  =  H*F(I) 

TEND(I)  =  T0(I)+TWRK(1,I)/4.0D0 
10  CONTINUE 
C 

c -  get  SECOND  ESTIMATE  - 

C 

CALL  DERIV6 (X0+H/4.0D0, TEND , F , HT, DELO , DELE , K , Kl , K2 ) 

DO  20  I  =  1,2 

TWRK(2,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(TWRK(1,I)*3,0D0+TWRK(2,I)*9.0D0)/32,0D0 
20  CONTINUE 
C 

c -  get  third  estimate  - 

c 

CALL  DERIV6 (X0+3.0D0*H/8. ODO, TEND, F,HT, DELO, DELE, K,K1,K2) 

DO  30  I  =  1,2 

TWRK(3,I)  =  H*F(I) 

TEND(I)  =  T0(I)+(TWRK(1,I)*1932.0D0-TWRK(2,I)*7200.0D0 
&  +  TWRK(3,I)*7296.0D0)/2197.0D0 

30  CONTINUE 
C 

c -  get  fourth  estimate  - 

c 

CALL  DERIV6{X0+12 .0*H/ 13.0, TEND, F,HT, DELO, DELE, K,K1,K2) 

DO  40  I  =  1,2 

TWRK(4,1)  =  H+F(I) 


139 


TEND(I)  =  T0(I)+(439.ODO*TVmK(l,I)/216.OD0-8.0D0*TWRK(2,I) 
&  +  3680.0DO*TWRK(3,I)/513.0DO-845.0DO*TWRK(4,I)/4104.0DO) 

40  CONTINUE 
C 

C -  GET  FIFTH  ESTIMATE  - 


C 

CALL  DERIV6 ( XO +H , TEND , F , HT , DELO , DELE , K , K 1 , K2 ) 

DO  50  I  =  1,2 

TWRK(5,I)  =  H*F(I) 

TEND(I)  =  TO(I)-8.0DO*TWRK(1,I)/27-0+2.0DO*TWRK(2,I) 

&  -  3544.0DO*TWRK(3,I)/2565.0DO+1859.0DO*TWRK(4,I)/4104.0DO 

&  -  11.0D0*TWRK(5,I)/40-0D0 

50  CONTINUE 
C 

C -  GET  SIXTH  ESTIMATE  - 

C 

CALL  DERIV6 ( XO+H/2 . ODO , TEND, F, HT, DELO , DELE , K, Kl , K2 ) 

DO  60  I  =  1,2 

TWRK(6,I)  =  H*F(1) 

60  CONTINUE 

C 

C -  ESTIMATE  ERROR  BY  COMPUTING  DIFFERENCE  BETWEEN  FOURTH  AND 

C  FIFTH  ORDER  EQUATIONS 

C 

SUM  =  O.ODO 
ERROR  =  O.ODO 
DO  70  I  =  1,2 

SUM  =  DABS<TWRK(1,I)/360.0D0-128.0D0*TWRK(3,I)/4275.0D0 

&  -  2197. 0D0*TWRK(4,I)/75240.0D0+TWRK(5,I)/50. ODO 

&  +  2.0D0*TWRK(6,I)/55.0D0) 

ERROR  =  DMAXl( ERROR, SUM) 

70  CONTINUE 

C 

c -  IF  ERROR  LESS  THAN  TOLERANCE,  COMPUTE  X  AT  END  OF  - 

C  INTERVAL  FROM  A  WEIGHTED  AVERAGE  OF  SIX  ESTIMATES 

C 

IF  ( ERROR. LT.TOL)  THEN 
DO  80  I  =  1,2 

TEND(l)  =  T0(1)+16,0D0*TWRK(1,I)/135.0D0 


&  +  6656. 0D0*TWRK(3,I)/12825. ODO 

&  +  28561. 0D0*TWRK(4,I)/56430. 0-9. 0D0*TWRK(5,I)/50. ODO 

&  +  2.0D0*TWRK(6,I)/55.0D0 

T0(I)  *  TEND(I) 

80  CONTINUE 

XO  =  XO+H 
END  IF 
C 

c -  IF  ERROR  GREATER  THAN  TOLERANCE,  REDUCE  STEP  AND  GO  - 

C  AGAIN 

C 
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IF ( ERROR. GT.TOL)  THEN 
H  =  H/2.0D0 
ENDIF 
C 

C -  IF  ERROR  IS  SIGNIFICANLTY  LESS  THAN  TOLERANCE,  RELAX  STEP 

C 

IF(ERROR.LT.H*TOL/10.0D0)  THEN 
H  =  H*2.0D0 
ENDIF 
C 

C -  IF  OVERSHOOT,  REDUCE  STEP  - 

C 

IF(XO+H.GT.XEND)  THEN 
H  =  XEND-XO 
ENDIF 
C 

ENDIF 

C 

C 

GO  TO  1 
END 
C 

C7E******  TRSY:  DERIV6 (X,T,F,HT,DEL0,DELE,K,K1,K2)  ****************** 
C 

C  PARAMETERS ; 

C 

C  X  -  INDEPENDENT  VARIABLE 

C  T  -  THE  ARRAY  THAT  HOLDS  THE  INITIAL  VALUES  OF  THE  FUNCTIONS 

C  F  -  AN  ARRAY  THAT  HOLDS  VALUES  OF  THE  DERIVATIVES 

C  HT  -  HEIGHT  OF  FIN 

C  DELO  -  FIN  THICKNESS  AT  BASE 

C  DELE  -  FIN  THICKNESS  AT  TIP 

C  K  -  FIN  CONDUCTIVITY 

C  K1  -  CONSTANT  =  2.0*SB*EMIS 

C  K2  -  CONSTANT  =  E*ABS 

C 

C - 

c 

SUBROUTINE  DERI V6 ( X , T , F , HT , DELO , DELE , K , K 1 , K2 ) 

REAL*8  X,T(2) ,F(2) ,HT, DELO, DELE, K,K1,K2, DEL, DELP 
DEL(X)  =  DELO+(DELE-DELO)*X/HT 
DELP  =  (DELE -DELO )/HT 
F(l)  =  T(2) 

F(2)  =  -(DELP*T(2) ) /DEL{X)+(K1*T(1)**4.0D0-K2) / (K*DEL(X) ) 

RETURN 

END 

C 

C7F  *****  TRSY:  MDLIN6 (FCN6,X1,X2,XR,XTOL,FTOL,NLIM,TB,TE,HT,DELO, 

C  DELE,K,K1,K2,Q,L) 

C 
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C  PARAMETERS : 

C 

C  FCN6  -  FUNCTION  THAT  COMPUTES  VALUES  FOR  F.  MUST  BE  DECLARED 
C  EXTERNAL  IN  CALLING  PROGRAM 

C  X1,X2  -  INITIAL  VALUES  OF  X.  F(X)  MUST  CHANGE  SIGNS  AT  THESE 
C  POINTS 

C  XR  -  RETURNS  THE  ROOT  TO  THE  MAIN  PROGRAM 

C  XTOL  -  TOLERANCE  FOR  X 

C  FTOL  -  TOLERANCE  FOR  F 

C  NLIM  -  LIMIT  TO  NUMBER  OF  ITERATIONS 

C  TB  -  BASE  TEMPERATURE 

C  TE  -  TIP  TEMPERATURE 

C  HT  -  FIN  HEIGHT 

C  DELO  -  BASE  THICKNESS 

C  DELE  -  TIP  THICKNESS 

C  K  -  FIN  CONDUCTIVITY 

C  K1  -  CONSTANT  =  2.0*SB*EMIS 

C  K2  -  CONSTANT  =  E*ABS 

CO-  HEAT  INPUT  THRU  BASE 

C  L  -  FIN  LENGTH 

C 

c - 

c 

SUBROUTINE  MDLIN6 ( FCN6 , XI , X2 , XR, TB, TE, HT, DELO , DELE, K, Kl , K2 , Q, L) 
REAL* 8  XERR, FSAVE , Fl , F2 , FR , Xl , X2 , XR, XTOL , FTOL , TB , TE , HT , DELO , 

&  DELE ,K,K1,K2,Q,L, FCN6 , XlSAV, X2SAV, FlSAV, F2SAV 

INTEGER  I, J, PASS 
CHARACTER* 2  ANS 

DATA  XTOL , FTOL/ 0 .0001D0,0.0000lD0/ 

C 

C -  INITIALIZE  VALUES  - 

C 

PRINT  200 
XlSAV  =  XI 
X2SAV  =  X2 

Fl  =  FCN6 (Xl,TB,TE,HT, DELO, DELE, K,K1,K2,Q,L) 

F2  =  FCN6(X2,TB,TE,HT, DELO, DELE, K,K1,K2,Q,L) 

FlSAV  =  Fl 
F2SAV  =  F2 
FSAVE  =  Fl 
C 

C -  INITIALIZE  COUNTER  - 

C 

1  =  0 
J  =  1 
PASS  =  0 
C 

C -  LIMIT  SEARCH  TO  INCREMENTS  OF  50  ITERATIONS  - 

C 


5  1  =  1  +  1 
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IF(I.GT. J*50)  GOTO  10 
XR  =  X2-F2*(X2-X1)/(F2-F1) 

FR  =  FCN6(XR,TB,TE,HT,DEL0,DELE,K,K1,K2,Q,L) 
IF(PASS.EQ.O)  PRINT  100,1 
IF(PASS.EQ.l)  PRINT  200 
IF(PASS.GE. 1)  PRINT  300,I,XR,FR 
IF(PASS.GE.  1)  PASS  =  PASS  1 
C 

C -  CHECK  STOPPING  CRITERIA  - 

C 

XERR  =  DABS(X1-X2)/2.0D0 

IF(XERR.LE.XTOL.OR.DABS(FR) .LE.FTOL)  RETURN 
C 

c - FIND  NEW  POINT - 

C 

IF(FR*F1.GE.0.0D0)  THEN 
XI  =  XR 
FI  =  FR 

IF(FR*FSAVE.GT.0.0D0)  F2  =  F2/2.0D0 
FSAVE  =  FR 
ELSE 
X2  =  XR 
F2  =  FR 

IF(FR*FSAVE.GT.0.0D0)  Fl  =  F1/2.0D0 
FSAVE  =  FR 
END  IF 
GOTO  5 
C 

c -  fail  to  find  root,  continue  search  7  - 

c 

10  PRINT  400,J*50 
READ  500,ANS 

IF(ANS.NE. * Y' .AND.ANS.NE. ’N' )  THEN 
PRINT  600 
GOTO  10 
ENDIF 

IF(ANS.EQ. 'N' )  GOTO  15 
J  =  J  +  1 
1  =  1-1 
GOTO  5 
C 

C -  FAIL  TO  FIND  ROOT,  LOOK  AT  FUNTION  VALUES  7 

C 

15  PRINT  700 
READ  500, ANS 

IF(ANS.NE. • Y' .AND.ANS.NE. 'N' )  THEN 
PRINT  600 
GOTO  15 
ENDIF 

IF(ANS.EQ. ’N’ )  STOP 
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XI  =  XISAV 
X2  =  X2SAV 
FI  =  FlSAV 
F2  =  F2SAV 
1  =  0 
J  =  1 
PASS  =  1 
GOTO  5 

C 

C  - 

C 


100  FORMAT('+',  ‘ITERATION  :  *,I3,‘  ') 

200  FORMAT (/,•  COMMENCING  SEARCH  FOR  ROOT  IN  INTERVAL _  ',/) 

300  FORMAT (•  AT  ITERATION *, 13 , 3X, ‘  X  =  •,D15.6,3X, ‘  F(X)  =  ‘,015.6) 

400  FORMAT (/,‘  FAILED  TO  FIND  ROOT  AFTER  ‘,13, ‘  ITERATIONS. ‘  , 

&  ‘  CONTINUE  SEARCH  (Y  OR  N)  7‘  ) 

500  FORMAT (A2) 

600  FORMAT(/‘  *****************  INCORRECT  ANSWER  *******************  ‘  ) 
700  FORMAT (/,•  WOULD  YOU  LIKE  TO  SEE  FITOCTION  VALUES  (Y  OR  N  )7  ‘,/, 

&  •  THIS  WILL  GIVE  SOME  IDEA  OF  WHAT  THE  FUNCTION  IS  DOING.  ‘,/, 

&  ‘  IF  F(XR)  IS  OSCILLATING  ABOUT  A  CERTAIN  VALUE  OR  ITS  ‘,/, 

&  •  ABSOLUTE  VALUE  IS  INCREASING  THEN  THERE  IS  LITTLE  • , / , 

&  ‘  PROBABILITY  THAT  A  ROOT  WILL  BE  FOUND.  ‘) 


END 

C 

C8A  *****  SI:  (TB, TE,L,HT, DEL, DELO, DELE, K, E, ,QI,Q)  ****************** 
C 

C  THIS  SUBROUTINE  CONVERTS  ARGUMENTS  TO  SI  UNITS 
C 

SUBROUTINE  S I ( TB , TE , L , HT , DEL , DELO , DELE , K , E , QI , Q ) 

REAL*8  X, TB,TE,L,HT, DEL, DELO, DELE, K,E,QI,Q 
REAL *8  DEGCF,MFT,WBTU,COND 


C 

C -  DEFINE  CONVERSION  FUNCTIONS  - 

C 

C -  CONVERT  TO  CENTIGRADE  FROM  FAHRENHEIT 

C 


DEGCF(X)  =  5.0D0*(X-32.0D0)/9.0D0 
C 

C - CONVERT  TO  METERS  FROM  FEET - 

C 

MFT(X)  =  X/3.280833D0 
C 

C - CONVERT  TO  W  FROM  BTU/HR - 

C 

WBTU(X)  =  X/3.412D0 
C 

C -  CONVERT  TO  W/M-K  FROM  BTU/HR-FT-F  (CONDUCTIVITY) 

C 


COND(X)  =  X*1.729577D0 
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c 

C -  CONVERT  ELEMENTS  - 

C 

TB  =  DEGCF(TB) 

TE  =  DEGCF(TE) 

L  =  MFT(L) 

HT  =  MFT(HT) 

DEL  =  MFT(DEL) 

DELO  =  MFT(DELO) 

DELE  =  MFT(DELE) 

K  =  COND(K) 

E  =  WBTU(E) 

QI  =  WBTU(QI) 

Q  =  WBTU(Q) 

RETURN 

END 

C 

C9A  *****  ENG:  (TB, TE,L,HT, DEL, DELO, DELE, K, E, ,QI,Q)  ****************** 
C 

C  THIS  SUBROUTINE  CONVERTS  ARGUMENTS  TO  ENG  UNITS 
C 

SUBROUTINE  ENG ( TB , TE , L , HT , DEL , DELO , DELE , K , E , QI , Q ) 

REAL*8  X, TB,TE,L,HT, DEL, DELO, DELE, K,E,QI,Q 
REAL*  8  DEGFC , FTM , BTUW , COND 


C 

C -  DEFINE  CONVERSION  FUNCTIONS  - 

C 

C -  CONVERT  TO  FAHRENHEIT  FROM  CENTIGRADE 

C 

DEGFC (X)  =  9.0D0*X/5.0D0  +  32.0D0 
C 

C - CONVERT  TO  FEET  FROM  METERS - 

C 

FTM(X)  =  X*3.280833D0 
C 

C - CONVERT  TO  BTU/HR  FROM  W - 

C 


BTUW(X)  =  X*3.412D0 


C 

C -  CONVERT  TO  BTU/HR-FT-F  FROM  W/M-K  (CONDUCTIVITY) 

C 

COND(X)  =  X/1.729577D0 
C 

C -  CONVERT  ELEMENTS  - 

C 

TB  =  DEGFC (TB) 

TE  =  DEGFC (TE) 

L  =  FTM(L) 

HT  =  FTM(HT) 


DEL  =  FTM (DEL) 
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DELO  =  FTM(DELO) 
DELE  =  FTH(DELE) 
K  =  COND(K) 

E  =  BTUW(E) 

QI  =  BTUW(Ql) 

Q  =  BTUW(Q) 

RETURN 

END 


********************************************************************** 

********************************************************************** 
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